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Abstract 

The application of nonlinear dynamics to improve the understanding of numerical uncer- 
tainties in computational fluid dynamics (CFD) is reviewed. Elementary examples in the use 
of dynamics to explain the nonlinear phenomena and spurious behavior that occur in numerics 
are given. The role of dynamics in the understanding of long time behavior of numerical 
integrations and the nonlinear stability, convergence, and reliability of using time-marching 
approaches for obtaining steady- state numerical solutions in CFD is explained. The study is 
complemented with examples of spurious behavior observed in CFD computations. 
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I. Introduction 


The authors’ view and experience in the application of nonlinear dynamics and bifurcation 
theory to improve the understanding of numerical uncertainties and their effects in computational 
fluid dynamics (CFD) are reviewed. The use of dynamics to illuminate how numerical methods 
work for strongly nonlinear problems is indirectly addressed. Simple nonlinear model equations 
are used to illustrate how the recent advances in nonlinear dynamical system theory can provide 
new insights and further the understanding of nonlinear effects on the asymptotic behavior 
of numerical algorithms commonly used in CFD. The discussion is complemented with CFD 
examples containing spurious behavior (numerical artifacts) in steady and unsteady flows. This 
topic is part of a long term research referred to as the “Dynamics of Numerics for CFD”. Here 
“dynamics” is used loosely to mean the dynamical behavior of nonlinear dynamical systems 
(continuum or discrete) and “numerics” is used loosely to mean the numerical methods and 
procedures in solving dynamical systems. For this paper, the phrase “to study the dynamics of 
numerics” (dynamical behavior of a numerical scheme) is restricted to the study of nonlinear 
and long time behaviors of nonlinear difference equations resulting from finite discretizations 
of a nonlinear differential equation (DE) subject to the variation of discretized parameters such 
as the time step, grid spacing, numerical dissipation coefficient, etc. This topic belongs to 
a subset of a rather new field in numerical analysis and dynamical system theory sometimes 
referred to as “The Dynamics of Numerics and The Numerics of Dynamics,” named after 
the First IMA Conference on Dynamics of Numerics and Numerics of Dynamics, University of 
Bristol, England, July 3 1 - August 2, 1990. We emphasize here that in the study of the dynamics 
of numerics, unless otherwise stated, we always assume the continuum (governing equations) 
is nonlinear. Although this paper is intended primarily for computational fluid dynamicists, it 
can be useful for computational scientists, physicists, engineers and computer scientists who 
have a need for reliable numerical simulation. 

Since the late 1980’s, many CFD related journals imposed an editorial policy statement 
on numerical uncertainty which pertained mainly to the accuracy issue. However, the study 
of numerical uncertainties in practical computational physics encompasses very broad subject 
areas. These include but are not limited to (a) problem formulation and modeling, (b) type, 
order of accuracy, nonlinear stability, and convergence of finite discretizations, (c) limits 
and barriers of existing finite discretizations for highly nonlinear stiff problems with source 
terms and forcing, and/or for wave propagation phenomena, (d) numerical boundary condition 
procedures, (e) finite representation of infinite domains (f) solution strategies in solving 
the nonlinear discretized equations, (g) procedures for obtaining the steady-state numerical 
solutions, (h) grid quality and grid adaptations, (i) multigrids, and (j) domain decomposition 
(zonal or multicomponent approach) in solving large problems. See the last six years of AIAA 
conference papers on numerical uncertainties in CFD and guidelines for code verification, 
validation and certification. See, for example, Mehta (1995), Melnik et al. (1994), Cosner et 
al. (1995), Demuren & Wilson (1994) Marvin (1993), Marvin & Holst (1990) and references 
cited therein. At present, some of the numerical uncertainties can be explained and minimized 
by traditional numerical analysis and standard CFD practices. Highly nonlinear and/or stiff 
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problems, however, do not always lend themselves to such treatment. At the same time, the 
understanding of the dynamics of numerics, in general, is at a very early stage of development 
and much remains to be learned. More theoretical development and extensive numerical 
experimentation are needed. Nevertheless, we believe that this approach can improve the 
understanding of numerical uncertainties and computational barriers in CFD in general, and 
in particular in combustion, direct numerical simulations, high speed and reacting flows, and 
certain turbulence models in compressible Navier-Stokes computations. 

A major stumbling block in genuinely nonlinear studies is that unlike the linear model 
equations used for conventional stability and accuracy considerations in time-dependent partial 
differential equations (PDEs), there is no equivalent unique nonlinear model equation for 
nonlinear hyperbolic and parabolic PDEs for fluid dynamics. A numerical method behaving 
in a certain way for a particular nonlinear DE (PDE or ordinary differential equation (ODE)) 
might exhibit a different behavior for a different nonlinear DE even though the DEs are of 
the same type. See Jackson (1989) for the definition of genuinely (or strongly) nonlinear 
problems. On the other hand, even for simple nonlinear model DEs with known solutions, 
the discretized counterparts can be extremely complex, depending on the numerical methods. 
Except in special cases, there is no general theory at the present time to characterize the various 
nonlinear behaviors of the underlying discretized counterparts. Most often, the only recourse 
is a numerical approach. Under this constraint, whenever analytical analysis of the discretized 
counterparts is not possible, the associated dynamics of numerics such as bifurcation phenomena 
and asymptotic behavior are obtained numerically using supercomputers. The term “discretized 
counterparts’’ is used to mean the finite difference equations (or “discrete maps”) resulting 
from finite discretizations of the underlying differential equations (DEs). Analysis here means 
analysis of the discretized counterparts by (a) available theory (b) searching for asymptotic 
solution behavior without resorting to purely numerical computations and (c) using continuation 
types of approaches to trace out the bifurcation diagrams (Keller 1977). This is the case for 
most of the illustrations throughout this paper. It is hoped that we can encourage numerical 
analysts to construct practical algorithms (to avoid spurious dynamics) based on the numerical 
phenomena observed using supercomputers to balance advances of computations and analyses. 
We also hope that it will strengthen the interface of numerical analysis with practical CFD 
applications and motivate CFD researchers who are looking for new. approaches and solutions 
to new or old but challenging problems. 

Due to the rapid independent development of dynamics of continuum and discrete maps, 
and the numerics of dynamics and the dynamics of numerics, it is extremely difficult to write 
a review paper on the subject. Topics discussed and references cited are only representative of 
the subject and reflect the authors’ experiences and preference toward certain areas at the time 
of this writing. Only selected properties of the observed phenomena are discussed. 

Outline: A rather detailed background, motivation and subtleties of the subject will be given 
in Section II due to the relatively new yet interdisciplinary nature of this research topic for 
CFD. The subtleties discussed are necessary to illuminate and isolate the sources of numerical 
uncertainties due to factors such as slow convergence or nonconvergence of numerical schemes. 
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and nonlinear behavior of high-resolution shock-capturing schemes. These include but are 
not limited to spurious nonlinear behavior of numerics such as “spurious chaos”, “spurious 
traveling waves”, “spurious chaotic transient” in transition to turbulence flows, “spurious 
steady-state numerical solutions” and “spurious asymptotes” (e.g., spurious limit cycles). 
Here “spurious numerical solutions (and asymptotes)” is used to mean numerical solutions 
(asymptotes) that are solutions (asymptotes) of the discretized counterparts but are not solutions 
(asymptotes) of the underlying DEs. Asymptotic solutions here include steady-state solutions, 
periodic solutions, limit cycles, chaos and strange attractors. See Thompson & Stewart (1986) 
and Hoppensteadt (1993) for the definition of chaos and strange attractors. The background 
material includes the connection between nonlinear dynamics and fluid dynamics, between 
nonlinear dynamics and CFD, between nonlinear dynamics and time-marching approaches. 

Section III reviews the basic terminology in nonlinear dynamics and reviews selected 
examples from our previous work. These examples consist of nonlinear model ODEs and PDEs. 
Particular attention is paid to the isolation of the different nonlinear behavior and spurious 
dynamics due to some of the numerical uncertainties mentioned earlier. The numerical schemes 
considered are selected to illustrate the following different spurious behavior of the dynamics 
of numerics. 

(a) Occurrence of stable and unstable spurious asymptotes above the linearized stability limit 
of the scheme (for constant step sizes) 

(b) Occurrence of stable and unstable spurious steady states below the linearized stability 
limit of the scheme (for constant step sizes) 

(c) Interplay of initial data and time steps on the permissibility of spurious asymptotes 

(d) Linearized behavior vs. nonlinear behavior of numerical solutions 

(e) Stabilization of unstable steady states by super-stable implicit methods 

(f) Interference with the dynamics of the underlying implicit scheme by procedures in solving 
the nonlinear algebraic equations (resulting from implicit discretization of the continuum) 

(g) Dynamics of the linearized implicit Euler scheme vs. Newton’s method 

(h) Local error control in ODE solvers and the conferrance of global properties of the 
nonlinear ODEs 

(i) Spurious chaos and chaotic transients 

(j) Spurious dynamics independently introduced by spatial and time discretizations 

Sections IV and V illustrate examples of CFD computations that exhibit spurious behavior 
due to numerics. The discussion is divided into transient and steady-state computations with 
several examples for each category. Sections 4.2.1, 5.1, 5.2 and 5.4 were written by the 
original contributors of the respective work. Section IV was taken from Yee & Sweby (1996). 
It is mainly concerned with convergence rate and spurious behavior of time-marching to the 
steady states of high-resolution shock-capturing methods. Section V is concerned with a 
“numerically induced chaotic transient” computation, and spurious behavior and convergence 
rate of transient computations of high-resolution shock-capturing schemes. The final section 
discusses how to use existing tools in bifurcation theory to avoid convergence to the wrong 
steady states or asymptotes. These tools are based on the combined knowledge of recent 
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advances in the dynamics of the physical equations and the dynamics of the underlying finite 
discretizations using existing tools in bifurcation theory and nonlinear dynamics. Note that the 
reference list contains more references than are cited inside the text. These extra references are 
intended for interested readers to pursue this subject further. 


II. Background & Motivations 

Starting in the late 1970’s, there have been important developments and breakthroughs 
concerned with the theory of nonlinear dynamical systems. There was also an explosion of 
journal and conference papers, texts and reference books on the subject in the 1980’s and early 
1 990’s. See for example, Guckenheimer & Holmes ( 1 983), Thompson & Stewart ( 1 988), Seydel 
(1988), Hales & Kocak (1991), Stewart (1990), Wiggins (1990), and Hoppensteadt (1993). 
During the early 1980’s, a new area of applied mathematics emerged from the interaction 
of dynamical system theory and numerical analysis. These developments addressed mainly 
mathematical principles and their applications of numerics in the understanding of the dynamics 
of DEs without discussing the connection between dynamics and numerics for initial value 
problems (IVPs) and initial boundary value problems (IBVPs). There was, however, some 
discussion on this connection for boundary value problems (BVPs) (Doedel & Beyn 1981, 
Doedel 1986, Shubin et al. 1981, Kellogg et al. 1980, Peitgen et al. 1981 and Shreiber & Keller 
1983). Studies of BVPs of the elliptic type continue to the present day. See, for example, the 
SIAM Conference on Dynamical Systems, October 15-19, 1992, and the Proceedings of IMA 
Conferences on Dynamics of Numerics and Numerics of Dynamics, July 31 - Aug. 2, 1990, 
Bristol, England. 

Why is it important to understand the connection between dynamics and numerics for BVPs, 
IVPs and IBVPs? It stems from the fact that it is a standard practice to use numerics to discover 
dynamical properties of continuous systems. As a matter of fact, much of what we know about 
specific dynamical systems is usually obtained from numerical experiments. One not only can 
visualize the dynamics and the bifurcation phenomena associated with numerics, but in most 
of the cases the dynamical behaviors of the DEs are not amenable otherwise. Consequently, 
developments concerned with the connection between dynamics and numerics are necessary 
to bridge the gap in a better understanding of the dynamics of numerics and the numerics of 
dynamics. 

In the late 1980’s, developments concerned with the connection between dynamics and 
numerics for IVPs and IBVPs slowly emerged. See for example, Mitchell & Griffiths (1985), 
Griffiths & Mitchell (1988), Pruffer (1985), Iserles (1988), Lorenz (1989) and Sanz-Sema 
(1985, 1990), and Salas et al. (1986). These developments raised many interesting and 
important issues of concern that are useful to practitioners in computational sciences. The 
following lists some of the issues: 

(a) Can recent advances in dynamical systems provide new insights into better understanding 
of numerical algorithms and the construction of new ones? 
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(b) Can these advances aid in the determination of a more reliable criterion on the use of 
existing numerical schemes for strongly nonlinear problems? 

(c) Under what conditions should a nonlinear problem be treated as a genuinely nonlinear 
problem rather than as a simplified linear problem? 

(d) Do traditional convergence and linear stability analyses apply to asymptotic nonlinear 
behavior and in particular to long time numerical simulation of nonlinear evolutionary PDEs? 

(e) Does error control suppress spuriosity? To what extent does local error control confer 
global properties in IVPs (long time integration) of nonlinear DEs? 

(f) What is the influence of finite time steps and grid spacings rather than time steps and 
grid spacings approaching zero on the overall nonlinear behavior and stability of the scheme in 
terms of allowable initial data and discretized parameters? 

(g) How different is the dynamical behavior of different procedures in solving the nonlinear 
algebraic equations resulting from using implicit time discretizations? 

Since the early 1990’s, the use of dynamics to address long time behavior of numerical 
schemes for IVPs and IB VPs began to flourish. The more recent work includes the Conference 
on Dynamics of Numerics and Numerics of Dynamics (University of Bristol, July 31 - 
August 2, 1990), the Chaotic Numerics Workshop (Deakin University, Geelong, Australia, 
July 12-16, 1993), the Conference on Dynamical Numerical Analysis (Georgia Institute of 
Technology, Atlanta, Georgia, December 14-16, 1995), and the “Innovative Time Integrators 
Workshop’ ’ (Center for Mathematics and computer Science, Amsterdam, November 6-8, 1996, 
the Netherlands). These conferences were devoted almost entirely to dynamical numerical 
analysis. See the proceedings and references cited therein. See also Stuart (1994, 1995), 
Humphries (1992), Hairer et al. (1989), Aves et al. (1995), Corless (1994a, b), Dieci & Estep 
(1991), and Poliashenko & Aidun (1995). The majority of the later developments concentrated 
on long time behavior of ODE solvers using variable step size based on local error controls 
(Butcher 1987). This type of local error controls enjoyed much success in controlling accuracy 
and stability for transient computations. 

On the other hand, even though standard practice in ODE solvers uses local error control 
for the selection of step size, practical CFD applications rarely employ this type of variable 
step size control due to the lack of efficient and practical theory for highly coupled nonlinear 
time-dependent PDEs. See Section 3.6 for a discussion. It should also be noted that finite 
element methods use ODEs solvers for the time integration part of the solution. There is also 
some preliminary work on computability and error control in finite element methods. See 
Johnson (1995) and Johnson et al. (1995) and references cited therein. Indirectly, this type of 
approach does make use of adaptive step size error controls. There still remains the question of 
spurious dynamics due to spatial discretizations. In practical CFD computations, one usually 
uses fixed grid spacings and time step constraints based on a linearized stability requirement or 
the Courant-Friedrich-Lewy (CFL) condition. Usually, after the initial transient dies down, the 
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step sizes are nearly constant from one step to the next in time marching to the steady state. 

The caveat is that regardless of whether finite difference (and finite volume) or finite 
element methods are employed, when time-marching approaches are used to obtain steady-state 
numerical solutions, local error controls similar to that used in ODE solvers that were designed 
for accuracy purposes are neither practical nor appropriate to use, since such local step size error 
control methods might prevent the solution from reaching the correct steady-state solutions 
within a reasonable number of iterations. It is remarked that the standard practice of using 
“local time step” (varied from grid point to grid point with the same CFL) in time-marching 
to the steady state is not the same as the variable step size based on local error controls. 

The authors believe that the understanding of the dynamics of numerics for fixed step size is 
necessary from that aspect. Besides, the study of the dynamics of ODE solvers using variable 
step size based on local error control requires the knowledge of the constant step size case 
(Aves et al. 1995). In a series of papers, Yee et al. (1991), Yee & Sweby (1994, 1995a,b), 
Sweby et al. (1990, 1995), Sweby & Yee (1992, 1994), and Lafon & Yee (1991, 1992) studied 
the dynamics of finite discretization for fixed (constant) time steps. The examples used in 
these papers were deliberately kept simple to permit explicit analysis. The approach was to 
take nonlinear model ODEs and PDEs with known explicit solutions (the most straight forward 
way of being sure what is ‘really’ happening), discretize them according to various standard 
numerical methods, and apply techniques from discrete dynamics to analyze the behavior of the 
discretized counterparts. Particular attention was paid to the isolation of the different nonlinear 
behavior and spurious dynamics due to some of the numerical uncertainties listed in Section 
I. These studies revealed much rich dynamical phenomena that we believe are useful for CFD 
and at the same time were not amenable with available theory at the time. The present paper 
alludes mainly to this aspect of the dynamics of numerics. 


2.1. Nonlinear Dynamics and Fluid Dynamics 

Most of the available fluid dynamics and CFD related texts and reference books describe the 
Euler and Navier-Stokes equations in differential form as coupled systems of nonlinear PDEs. 
These equations are rarely classified as dynamical systems. However, fluid dynamicists are 
often interested in how the flow behaves as a function of one or more physical parameters. Of 
particular interest to fluid dynamicists is locating the critical value of the physical parameter 
where the fluid undergoes drastic changes in the flow behavior. Some examples are the 
prediction of transition to turbulence or laminar instability as a function of the Reynolds 
number, flow separation and stall as a function of Reynolds number and angle of attack, 
rotorcraft vibration as a function of rotation speed and flight speed, the occurrence of shock 
waves as a function of the body shape and/or Mach number, and the formation of vortices, 
flutter, and other flow phenomena as a function of the angle of attack or other physical 
parameters. Another application is in the area of aiding the understanding of the topology of 
flow patterns (flow visualizations) of laboratory experiments, observable physical phenomena 
and numerical data. An additional important topic for CFD is the control and optimization 
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of dynamical systems. This involves the application of optimization and control theory to 
dynamical systems. Researchers are beginning to use these interdisciplinary ideas to study, for 
example, the control of turbulence, the control of vortex generation and/or shock waves, the 
control of vibration in rotorcraft, and the control of aerodynamic noise such as sonic boom and 
jet noise. 

The application of dynamical system theory to the study of spatio-temporal instabilities of 
aerodynamic and hydrodynamic flows and chaotic systems in fluid dynamics was discussed 
respectively in the 1994 and 1996 von Karman Institute for Fluid Dynamics Lecture Series. 
How the solution behaves as one or more of the system parameters is varied is precisely the 
definition of dynamical systems and bifurcation theory. According to Ian Stewart (1990) 

“Bifurcation theory is a method for finding interesting solutions to nonlinear equations by 
tracking dull ones and waiting for them to lose stability .” 

As evident from the Third SIAM Conference on the Application of Dynamical Systems, 
May 21-24, 1995, Snowbird, Utah, presentations in treating the various fluid flow equations as 
dynamical systems have pushed these topics to the forefront of applied mathematical research. 


2.2. Nonlinear Dynamics and CFD 


When we try to use numerical methods to gain insight into the fluid physics, there is an 
added new dimension to the overall problem. Even though we freeze the physical parameters 
of the governing equations, the resulting discretized counterparts (from finite discretizations of 
the governing equations) are not just a nonlinear system of difference equations, but are also 
a nonlinear but discrete dynamical system on their own. From nonlinear dynamics, we know 
that discrete dynamical systems possess much richer dynamical behavior than the continuum 
dynamical systems. These resulting discrete dynamical systems are a function of all of the 
discretized parameters which are not present in the governing equations. See Section III for a 
discussion. This is one of the key factors in influencing the numerical solution to depart from 
the physical ones if the governing equations are strongly nonlinear and stiff. 

Of course, before analyzing the dynamics of numerics, it is necessary to analyze (or 
understand) as much as possible the dynamical behavior of the governing equations and/or 
the physical problems using theories of ODEs, PDEs, dynamical systems of ODEs and PDEs, 
and also physical guidelines. In fact, a knowledge of the theories of ODEs, PDEs, dynamical 
behavior of nonlinear DEs, and the dynamical behavior of nonlinear discrete maps (difference 
equations) is a prerequisite to the study of dynamics of numerics. In an idealized situation, if 
one knows the dynamical behavior of the governing equations, one can then construct suitable 
numerical methods for that class of dynamical systems. Consequently, spurious dynamics due 
to numerics can be minimized and that computation and analysis kept pretty much in tune. 
However, as applied scientists want to push the envelope of understanding of realistic flows and 
configurations further, dependence on the numerics takes over even though rigorous analysis 
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lags behind. Starting in the late 1970’s the advances in computer power resulted in attempts 
to use CFD to replace wind tunnel experiments and use numerics to understand dynamical 
systems. The gap between computation and analysis increased. The nonlinear behavior of 
commonly used algorithms in CFD was not well understood but at the same time applied CFD 
increased the intensity of using these algorithms to solve more complex practical problems 
where the flow physics and configurations under consideration were not understood, and were 
either too costly for or not amenable to laboratory experiments. CFD was and remains in 
a stage where computation is ahead of analysis. In other words, we usually do not know 
enough about the solution behavior of the underlying DEs in practice, and we are at the stage 
where the understanding of the dynamics of the DEs and the understanding of the dynamics of 
numerics are in tandem, and they both are rapidly growing research areas. With this in mind, 
we summarize some of the sources of nonlinearities in the study of dynamics of numerics for 
CFD. 

The sources of nonlinearities that are well known in CFD are due to the physics. Examples 
of nonlinearities due to the physics are convection, diffusion, forcing, turbulence source terms, 
reacting flows, combustion related problems, or any combination of the above. The less familiar 
sources of nonlinearities are due to the numerics. There are generally three major sources: 

(a) Nonlinearities due to time discretizations — the discretized counterpart is nonlinear in 
the time step. Examples of this type are Runge-Kutta methods. It is noted that linear multistep 
methods (LMMs) are linear in the time step. See Lambert ( 1 973) for the forms of these methods. 

(b) Nonlinearities due to spatial discretizations — in this case, the discretized counterpart can 
be nonlinear in the grid spacing and/or the scheme. Examples of nonlinear schemes are the total 
variation diminishing (TVD) and essentially nonoscillatory (ENO) schemes. See Yee (1989) 
and references cited therein for the forms of these schemes. 

(c) Non linearities due to complex geometries, boundary interfaces, grid generation, grid 
refinements and grid adaptations — each of these procedures can introduce nonlinearities. 

The behavior of the above nonlinearities due to the numerics are not well understood. Only 
some preliminary development is beginning to emerge recently. 


2.3. Dynamics of Numerical Approximations of ODEs vs. Time- Dependent PDEs 

Recent analyses and studies have shown that spurious numerical solutions can be inde- 
pendently introduced by time and spatial discretizations. Take the case when the ODEs are 
obtained from semi-discrete approximations of PDEs, the resulting system of ODEs contains 
more parameters (due to spatial discretizations) as opposed to physical problems governed 
by ODEs. The parameters due to spatial discretizations for the semi-discrete approximation 
becomes the system parameter (instead of the discretized parameter) of the resulting system of 
ODEs. Depending on the differencing scheme the resulting discretized counterparts of a PDE 
can be nonlinear in A t, the grid spacing A* and the numerical dissipation parameters, even 
though the PDEs have only one parameter or none. One major consideration is that one might 
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be able to choose a “safe” numerical method to solve the resulting system of ODEs to avoid 
spurious stable numerical solutions due to time discretizations. However, spurious numerical 
solutions and especially spatially varying steady states introduced by spatial discretizations in 
nonlinear hyperbolic and parabolic PDEs for CFD applications appear to be more difficult to 
avoid due to the use of a fixed mesh. In the case of the semi-discrete approach such as methods 
of lines or finite element methods, if spurious numerical solutions due to spatial discretizations 
exist, the resulting ODE system has already inherited this spurious feature as part of the exact 
solution of the semi-discrete case. Thus care must be taken in using the ODE solver computer 
packages for PDE applications. See Lafon & Yee (1991, 1992) and Section 3.7 for a discussion. 


2.4. Nonlinear Dynamics and Time-Marching Approaches 


The use of time-marching approaches to obtain steady-state numerical solutions has been 
considered the method of choice in CFD for nearly two decades since the pioneering work of 
Crocco (1965) and Moretti & Abbett (1966). Moretti and Abbett used this approach to solve 
the inviscid supersonic flow over a blunt body without resorting to solving the steady form of 
PDEs of the mixed type. The introduction of efficient CFD algorithms of MacCormack (1969), 
Beam & Warming (1978), Briley & McDonald (1977) and Steger (1978) marked the beginning 
of numerical simulations of 2-D and 3-D Navier-Stokes equations for complex configurations. 
It enjoyed much success in computing a variety of weakly and moderately nonlinear fluid 
flow problems. For strongly nonlinear problems, the situation is more complicated. To aid 
the understanding of the scope of the situation, first, we have to identify all the sources of 
nonlinearities. Second, we have to isolate all elements and issues of numerical uncer tain ties 
due to these nonlinearities in time-marching to the steady state. The following isolates some of 
the key elements and issues of numerical uncertainties in time-marching to the steady state. 

Solving an IBVP with Unknown Initial Data '. In time-marching approaches, one transforms 
a BVP to an IBVP with unknown initial data. Elements such as reliability, convergence with 
unknown initial data and rapid convergence to the correct steady state are the predominant 
requirements. The time differencing in this case acts as a pseudo time. Although it is well 
known from linearized stability analysis that only a subset of the numerical solutions for 
certain ranges of the discretized parameters (these ranges can be continuous or disjoint for 
problems with source terms) and boundary conditions mimic the true solution behavior of the 
governing equations, it is less well known that outside these safe regions the numerical solution, 
depending on the initial data, do not necessarily undergo instabilities. In addition, there exist 
asymptotic numerical solutions that are not solutions of the continuum even inside the safe 
regions. Unlike nonlinear problems, the numerical solutions of linear or nearly linear problems 
are “independent” of the discretized parameters and initial data as long as the discretized 
parameters are inside the stability limit (CFL condition). We put “independent” in quotations 
here to mean that the numerical solutions behave the same up to the order of accuracy and 
grid spacing of the scheme. That is, the topological shapes of these solutions remain the same 
within the stability limit and accuracy of the scheme for linear behavior. Section 3.4 illustrates 
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the strong dependence of the numerical solution on initial data for nonlinear DEs. It turns out 
that if constant step sizes are used, stability, convergence rate and permissibility of spurious 
numerical solutions are intimately related to the choice of initial data (or start up solution). 

Reliability of Residual Test: The deficiency of the use of residual tests in detecting the 
convergence rate and the convergence to the correct steady-state numerical solutions is now 
briefly discussed. Consider a quasilinear PDE of the form 

u t = G(u, u. e), (2.1) 

where G is nonlinear in u, u, and u B ». The values a and e are system parameters. For 
simplicity, consider a two time level and a (p + q + 1 ) point grid stencil numerical scheme of 
the form 


u 


"+ 1 — „ n 


= <*” - A4,Ai) 


( 2 . 2 ) 


for the PDE (2.1). Note that the discussion need not be restricted to explicit methods or two 
time level schemes. Let U*, a vector representing (u^ +9 , ..., uj , ..., be a steady-state 

numerical solution of (2.2). When a time-marching approach such as (2.2) is used to solve 
the steady-state equation G(u,u m ,u mm ,a,e) = 0, the iteration typically is stopped when the 
residua] H and/or some li norm of the dependent variable « between two successive iterates is 
less than a pre-selected level. 


Aside from the various standard numerical errors such as truncation error, machine round-off 
error, etc., there is a more fundamental question of the validity of the residual test and/or £3 norm 
test. If the spatial discretization happens to produce spurious steady-state numerical solutions, 
these spurious solutions would still satisfy the residual and li norm tests in a deceptively smooth 
manner. Moreover, aside from the spurious solution issue and depending on the combination of 
time as well as spatial discretizations, it is not easy to check whether G(ti% u* , u*„, a, e) — » 0 
even though H(U*,a,e,At,Ax) — ► 0, since spurious steady states (and asymptotes) can be 
independently introduced by spatial and time discretizations. This is contrary to the ODE case, 
where if u* is a spurious steady state of du/dt = S(u), then S(u*) 7 ^ 0. If a steady state has 
been reached with a rapid convergence rate, it does not imply that the obtained steady state is 
not spurious. 


Methods Used to Accelerate Convergence Process: Methods such as iterations and relaxation 
procedures, and/or convergence acceleration methods such as conjugate gradient methods have 
been utilized to speed up the convergence process (Saad 1994). Also techniques such 
as preconditioning (Turkel 1993) and multigrid (Wesseling 1992) combined with iteration, 
relaxation and convergence acceleration procedures are commonly used in CFD. Depending on 
the type of PDEs, proper preconditioners can be established for the PDEs or for the particular 
discretized counterparts. Multigrid methods can be applied to the steady PDEs or the time- 
dependent PDEs. In either case, a combination of these methods can still be viewed as pseudo 
time-marching methods (but not necessarily of the original PDE that was under consideration). 
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However, if one is not careful, numerical solutions other than the desired one can be obtained 
in addition to spurious asymptotes due to the numerics. From here on the term “time-marching 
approaches” is used loosely to include all of the above. It is remarked that multigrid methods 
can be viewed as the (generalized) spatial counterpart of variable time step control in time 
discretizations. 

Method* in Solving the Nonlinear Algebraic Equations From Implicit Method < : When 

implicit time discretizations are used, one has to deal with solving systems of nonlinear 
algebraic equations. Aside from the effect of the different methods to accelerate the conver- 
gence process discussed previously, we need to know how different the dynamical behavior is 
for the different procedures (e.g., iterative vs. non-iterative) in solving the resulting nonlinear 
difference equations. See Yee & Sweby (1994, 1995a,b) and the next section for a discussion. 

Mismatch in Implicit Scheme s \ It is standard practice in CFD to use a simplified implicit 
operator (or mismatched implicit operators) to reduce CPUs and to increase efficiency. These 
mismatched implicit schemes usually consist of the same explicit operator but different 
simplified implicit operators. The implicit time integrator is usually of the LMM type. One 
popular form of the the implicit operator is the so called ‘ ‘delta formulation” (Beam & Wanning 
1978). The advantage of the delta formulation is that if a steady state has been reached, the 
numerical solution is independent of the time step. The original logic in constructing this 
type of scheme is that the implicit operators act as a relaxation mechanism. However, from a 
dynamical system standpoint, before a steady state is reached, the nonlinear difference equations 
representing each of these simplified implicit operators are different from each other. They have 
their own dynamics as a function of the time step, grid spacing and initial data. They also can 
exhibit different types of nonlinear behavior if one is solving strongly nonlinear time-dependent 
PDEs. Thus, even when steady states have been reached, they might not converge to the same 
steady state due to the different forms of the implicit operator. In addition, even if they converge 
to the same steady state, that steady state can still be spurious since the same explicit operator 
is used. This might be due to the fact that spatial discretizations can introduce spurious steady 
states. Consequently, these mismatched implicit operators can have different spurious dynamics 
and/or different convergence rates for the entire solution procedure. Section 4.3 describes some 
examples. 

Nonlinear Schemes; . It is well known that all of the TVD, total variation bounded (TVB) 
and ENO schemes (see Yee (1989) or references cited therein) are nonlinear schemes in the 
sense that the final algorithm is nonlinear even for the constant-coefficient linear PDE. These 
types of schemes are known to have a slower convergence rate than classical shock-capturing 
methods and can occasionally produce unphysical solutions for certain combinations of entropy 
satisfying parameters and flux limiters (in spite of the fact that entropy satisfying TVD, TVB 
and ENO schemes can suppress unphysical solutions). See Yee (1989) for a summary of the 
subject. The second aspect of these nonlinear schemes is that even if the numerical method is 
formally of more than first-order and if the approximation converges, the rate may still be only 
first-order behind the shock (not just around the shock). This can happen for systems where one 
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characteristic may propagate part of the error at a shock into the smooth domain. Engquist & 
Sjogreen (1996) illustrate this phenomena with examples. See Section 4.2 for a discussion. The 
third aspect of these higher-order nonlinear schemes is their true accuracy away from shocks. 
See Donat (1994), Casper & Carpenter (1995) and Section 5.4 for a discussion. 

Schemes That are Linear vs. Nonlinear in At . The obvious classification of time-accurate 
schemes for time-marching approaches to the steady state are explicit, implicit, and hybrid 
explicit and implicit methods. Although many of the added numerical procedures, discussed 
in the previous paragraph, used to help speed up convergence to steady states can apply to 
both explicit and implicit methods, there are distinct differences between the two. Usually, 
one tries to enlarge what is the equivalent of the linearized time step constraints imposed by 
the explicit schemes for rapid convergence. On the other hand, if implicit methods are used, 
unconditionally stable implicit methods are usually employed. Linearized stability constraints 
are not the problem. Efficient methods for solving the nonlinear algebraic equations which 
guarantee rapid convergence to the correct steady states are then the main concern. 

A less commonly known classification of numerical schemes for time-marching approaches 
is the identification of schemes that are linear or nonlinear in the time step (A t) parameter 
space when applied to nonlinear DEs. As mentioned before, all LMMs (explicit or implicit) are 
linear in A t and all multistage Runge-Kutta methods are nonlinear in At. Lax-Wendroff and 
MacCormack type of non-separable full discretizations also are nonlinear in At. A desirable 
property for a scheme that is linear in At is that, if the numerical solution converges, its steady- 
state numerical solutions are independent of the time step. On the other hand, the accuracy of 
the steady-state numerical solutions (also for time-accurate numerical solutions) depends on At 
if the scheme is nonlinear in At. Certain of these types of schemes are more sensitive to At 
than others. For example, Lax-Wendroff and MacCormack methods (MacCormack 1969) are 
more sensitive than the Lerat variant (Lerat & Sides 1988). A less known property of schemes 
that are nonlinear in At is that this type of scheme has an important bearing on the existence 
of spurious steady-state numerical solutions due to time discretizations. Although schemes like 
LMMs are immune from exhibiting spurious steady-state numerical solutions, as seen in Yee 
& Sweby (1994, 1995a,b), a wealth of surprisingly nonlinear behavior of implicit LMMs that 
had not been observed before were uncovered by the dynamical approach. See the next section 
for a review. 

Adaptive Time Step Based on Local Error Control . It is a standard practice in CFD to use 
“local time step” (varied from grid point to grid point using the same CFL) for nonuniform 
grids. However, except in finite element methods, adaptive time step based on local error 
control is rarely use in CFD. Adaptive time step is built in for standard ODE solver computer 
packages (Butcher 1987). It enjoyed much success in controlling accuracy and stability for 
transient (time-accurate) computations. The issue is to what extent does this adaptive local 
error control confer global properties in long time integration of time-dependent PDEs? Can 
one construct similar error control that has guaranteed and rapid convergence to the correct 
steady-state numerical solutions in the time-marching approaches for time-dependent PDEs? 
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Nonunique Steady- State Solutions of Nonlinear DEs v*. Spurious Asymptotes '. The phe- 
nomenon of generating spurious steady-state numerical solutions (or other spurious asymptotes) 
by certain numerical schemes is often confused with the nonuniqueness (or multiple steady 
states) of the governing equation. In fact, the existence of nonunique steady-state solutions of 
the continuum can complicate the numerics tremendously (e.g., the basins of attraction - which 
initial data lead to which asymptote) and is independent of the occurrence of spurious asymp- 
totes of the associated scheme. But, of course, a solid background in the theory of nonlinear 
ODEs and PDEs and their dynamical behavior is a prerequisite in the study of the dynamics 
of numerical methods for nonlinear PDEs. A full understanding of the subject can shed some 
light on the controversy about the “true” existence of multiple steady- state solutions through 
numerical experiments for certain flow types of the Euler and/or Navier-Stokes equations. 


III. Elementary Examples 

This section reviews the fundamentals of and illustrates with elementary examples, the 
spurious behavior of commonly used time discretizations in CFD. Except for Section 3.6, 
details of these examples can be found in our earlier papers. 

3.1. Preliminaries 

Consider an autonomous nonlinear ODE of the form 

^ = aS(u), (3.1) 

where a is a parameter and S(u) is nonlinear in u. Autonomous here means that i does not 
appear explicitly in the function S. For simplicity of discussion, we consider only autonomous 
ODEs where a is linear in (3.1); i.e., a does not appear explicitly in 5. 

A fixed point u* of an autonomous system (3.1) is a constant solution of (3.1); that is 

5(u*) = 0. (3.2) 

We remark that the terms “equilibrium points”, “critical points”, “singular points”, “sta- 
tionary points”, “asymptotic solutions” (we are excluding periodic solutions for the current 
definition), “steady-state solutions” and “fixed points” are sometimes used with slightly 
different meanings in the literature, for example, in bifurcation theory. However, for the current 
discussion and for the majority of the nonlinear dynamics literature, these terms are used 
interchangeably. Note that certain researchers reserve the term “fixed point” for discrete maps 
(difference equations) only. 

Consider a nonlinear discrete map from the finite discretization of (3.1) 
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u n+1 = u n + D(u",r), 


(3.3) 

where r = a At and D(u n ,r) is linear or nonlinear in r depending on the numerical scheme. 
Here the analysis is similar if (3.3) involves more than two time levels. Examples to illustrate 
the dependence on the numerical schemes for cases where D is linear or nonlinear in the 
parameter space will be given in a subsequent section. 

A fixed point u* of (3.3) (or fixed point of period 1) is defined by u" +1 = u", or 

u* = «* +D(u*,r), (3.4a) 

i.e., 

D(u*,r) = 0. (3.4b) 

One can also define a fixed point of period p (or periodic solution of period p), where p is a 
positive integer by requiring that u n+p = u n or 

u* = E p (u\r) but u*^V. r ) for 0 < k < p. (3.5) 

Here, E p (u*,r) means that we apply the difference operator E p times, where E(u n ,r) = 
u" -|- D(u n , r ). For example, a fixed point of period 2 means u n+3 = u n or 

u* = E(E(u\r)). (3.6) 

In the context of discrete systems, the term “fixed point” without indicating the period 
means “fixed point of period 1 ” or the steady-state solution of (3.3). Interchangeably, we also 
use the term asymptote to mean a fixed point of any period. 

In order to illustrate the basic idea, the simplest form of the Ricatti ODE, i.e., the logistic 
ODE with 


— = aS(u) = au(l-u) 


is considered. For this ODE, the exact solution is 


(3.7a) 


u(<) = 


u 


o 


u° + (1 - u°)e -at ’ 


(3.7b) 


where u° is the initial condition. The fixed points of the logistic equation are roots of 
u*(l — u*) = 0; it has two fixed points u* = 1 and u* = 0. 


To study the stability of these fixed points, we perturb the fixed point with a disturbance 
and obtain the perturbed equation 


f = qS ("' + () - 


(3.8) 
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Next, S(u* + £) can be expanded in a Taylor series around u*, so that 


dt 


= Q 


S(u-) + 54u-K + -S uu (a > )( 1 + 


(3.9) 


where 5 u (u*) = ^| u . . Stability can be detected by examining a small neighborhood of the 
fixed point provided that, for a given a, u* is a hyperbolic point (Seydel 1988) (i.e., if the 
real part of aS u (u*) / 0). Under this condition £ can be assumed small, its successive powers 
£ 2 , £ s , ... can normally be neglected and the following linear perturbed equation is obtained 


^ = aS u (u*)£. (3.10) 

The fixed point u* is asymptotically stable if a5«(u*) < 0 whereas u* is unstable if 
a5 u (u*) > 0. If a5 u (u*) = 0, a higher order perturbation is necessary. If u* is not a 
hyperbolic point, the behavior of (3.10) does not infer the behavior of the original unperturbed 
equation. 


For the logistic ODE, the fixed points are hyperbolic. Thus the linearized analysis suffices 
(i.e., the original equation has the same local behavior as the perturbed equation (3. 10)). If we 
perturb the logistic equation around the fixed point with a > 0, we find that u* = 1 is stable 
and u* = 0 is unstable. It is well known that the global asymptotic solution behavior of the 
logistic ODE is that for any u° > 0, the solution will eventually tend to u* = 1. Figure 3.1 
shows the asymptotic solution behavior of the logistic ODE. 

Now, let us look at three of the well known schemes for IVPs of ODEs. These are 
explicit Euler (Euler or forward Euler), leapfrog and Adams-Bashforth. For the ODE (3.1), the 
dynamical behavior of their corresponding discrete maps is well established. The explicit Euler 
method is given by 


u" +1 = u" + rS(u n ), (3.11) 

and after a linear transformation it is the well known logistic map (Hoppensteadt 1993). The 
leapfrog method can be written as 


u n+i = u n-i + 2r S(u n ), (3.12) 

and it is a form of the Henon map (Devaney 1987). The Adams-Bashforth method yields 


u 


n+1 ..n 


= U" + 


3S(u n ) - 5(u n_1 ) 


(3.13) 


again a variant of the Henon map that has been discussed by Pruffer (1985) in detail. 


We can determine fixed points of the discrete maps (3.11 )-(3. 1 3) and their stability properties 
in a manner similar to that for the ODE. It turns out that all three of the discrete maps have the 
same fixed points as the ODE (3.1) — a desired property which is important for obtaining the 
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correct steady states of nonlinear DEs numerically. An examination of (3.1 1 )-(3. 1 3) reveals 
that the discretized parameter r appears linearly in these discrete maps. As will be seen in the 
next section, a necessary condition for the occurrence of spurious steady states by any scheme 
is that the discretized parameter should appear nonlinearly in the underlying discrete map. 
Consequently the existence of spurious steady states is not possible for (3.1 1)-(3.13). 

The corresponding linear perturbed equation for the discrete map (3.3), found by substituting 
«" = u* + £ n in (3.3) and ignoring terms higher than £ n , is 

t n+1 = C [1 + A<D u (u*, A<)] . (3.14) 

Here the parameter a of the ODE has been absorbed in the parameter At based on the 
assumption that a does not appear explicitly in S(u). Depending on the scheme, D(u n , At) 
might be nonlinear in At. If nonlinearity in the parameter space At is introduced into the 
discretized counterpart, it increases the possibility that the dynamics of numerics deviates from 
the dynamics of the continuum. 

For stability we require 


|1 + AtI? u (ti*,At)| < 1. (3.15) 

A gain , for |1 + AtZ> u (u*, A<)| = 1, a higher order perturbation is necessary. For a fixed point 
of period p the corresponding linear perturbed equation and stability criterion are 


(" +p = C-E£K,A(), 

(3.16a) 

|£J(u’,A<)l <1, 

(3.16b) 

£'(<*”, Af) = £i!(u" +! '- , ,At)...£fi(o”,A()- 

(3.16c) 


For S(u) = u(l — u), the stability of the stable fixed points of period 1 and 2 for discrete 
maps (3.1 1)-(3.13) with r = a At are 

Explicit Euler: 

u* = 1 stable if 0 < r < 2 

period 2 stable if 2 < r < v/6 

Leapfrog: 

u* = 1 unstable for all r > 0 

chaotic solutions exist for all r no matter how small 
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Adams-Bashforth: 

u* = 1 stable if 0 < r < 1 

period 2 stable if 1 < r < \/2. 

Figure 3.2a shows the stable fixed point diagram of period 1,2, 4, 8 obtained for the 
explicit Euler scheme by solving numerically the roots of (3.1 1) (by setting u n+1 = u n ) for 
5(u) = tx(l - u). The r axis is divided into 1000 equal intervals. The numeric labeling of the 
branches denotes their period. The subscript “E” on the period 1 branch indicates the stable 
fixed point of the DE. From here on, r = aAt where a = a in all of the plots. The change of 
notation inside the plots is due to the plotting package. 

All of these three examples share a common property of not exhibiting spurious steady 
states. Two of these three examples serve to illustrate that operating with a time step beyond the 
linearized stability limit of the stable fixed points of the nonlinear ODEs does not always result 
in a divergent solution; spurious asymptotes of higher period can occur. This is in contrast to 
the ODE solution, where only a single stable asymptotic value u* = 1 exists for any a > 0 and 
any initial data u° > 0. The spurious asymptotes, regardless of the period, stable or unstable, 
are solutions in their own right of the discrete maps resulting from a finite discretization of the 
ODE. 


3.2. Spurious Asymptotic Numerical Solutions for Constant Step Sizes 

For the previous section we purposely picked the type of schemes that do not exhibit spurious 
fixed points, but allow spurious fixed points of period higher than 1. In this section numerical 
methods are purposely chosen so that the discretized parameter appears nonlinearly in the 
underlying discrete maps. Consequently, existence of spurious steady-state numerical solutions 
in these examples is possible. 


3.2.1. Explicit Methods 

Consider two second- and third-order Runge-Kutta schemes, namely, the modified Euler 
(R-K 2) and the improved Euler (R-K 2), Heun (R-K 3), Kutta (R-K 3), the fourth-order 
Runge-Kutta method (R-K 4), and two predictor-corrector methods (P-C 2 and P-C 3) (Lambert 
1973) of the forms 

Modified Euler (R-K 2) method: 

u n+1 = u n + rS fu n + , S n = S(u n ), (3.17) 
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Improved Euler (R-K 2) method: 


Heun (R-K 3) method: 


Kutta (R-K 3) method: 


R-K 4 method: 


u n+1 = u n + 


S n + S(u n + rS n )\, 



u n+1 =u n +^(ki+4k 2 +l fe,) 


kt = S " 




(3.18) 


(3.19) 


(3.20) 


(3.21) 
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Predictor-corrector for m = 2, 3 (P-C 2 & P-C 3): 


U< 0) = u n +rS r 


u 


(*+i) _ 


= u " + r 


u 


n+l _ ,,n 


= “ +2 


S n + 5 (fc) 

5" + 5 (m_1) 


fe = 0, 1, m — 1 


(3.22) 


3.2.2. Fixed Point Diagrams 

Using the same procedures as before, one can obtain the fixed points for each of the above 
schemes (3.17) - (3.22). Figures 3.2b - 3.2f show the stable fixed point diagrams of period 
1,2,4 and 8 for selected schemes for 5(u) = u(l - u). The unstable fixed points of any period 
are not plotted. See Yee et al. (1991) for the unstable fixed point diagrams. Some of the fixed 
points of lower period were obtained by closed form analytic solution and/or by a symbolic 
manipulator such as MAPLE (1988) to check against the computed fixed point The majority 
were computed numerically. The stability of these fixed points was examined by checking the 
discretized form of the appropriate stability conditions. The domain is chosen so that it covers 
the most interesting part of the scheme and ODE combinations, and is divided into 1000 equal 
intervals. In other words, spurious asymptotes may occur in other parts of the domain as well. 
The numeric labeling of the branches denotes their period, although some labels for period 4 
and 8 are omitted due to the size of the labeling areas. Again, the subscript “E” on the main 
period one branch indicates the stable fixed point of the DE while the subscript “S” indicates 
the spurious stable fixed points introduced by the numerical scheme. Spurious fixed points of 
period higher than one are obvious (since the ODEs under discussion only possess steady-state 
solutions) and are not labeled with a subscript “S”. Note that these diagrams, which for the 
most part appear to consist of solid lines, actually consist of points, which are only apparent 
in areas with high gradients. 

To contrast the results, similar stable fixed point diagrams were also computed for the ODE 

du 

— = qu(1-u)(6-u), 0 < & < 1, (3.23) 

that is, for a cubic nonlinearity for $(u) = u(l - u)(& - u). The stable fixed point for the ODE 
(3.23) in this case is u* = b and the unstable ones are u* = 0 and u* = 1. For any 0 < u° < 1 
and any a > 0, the solution will asymptotically approach the only stable asymptote of the ODE 
u* = b. 

By looking at the roots of the underlying discrete maps, it is readily realized that r appears 
nonlinearly in these discrete maps. In fact, the maximum number of stable and unstable fixed 
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points (real and complex) for each of the studied schemes (3.17)-(3.22) varied from 4 to 16 for 
5(u) = u(l — u) and 9 to 81 for S(u) = u(l — u)(b — u), depending on the numerical method 
and the r value. For certain r values, aside from the real steady states, these roots might be 
unstable and/or complex, but not for others. Fig. 3.3 shows the stable fixed point diagram by 
the modified Euler for four different values of b. 

Aside from the striking difference in topography in the stable fixed point diagrams of the 
above methods and ODE combinations, all of these diagrams (except the P-C 2 method) have 
one common feature: they all exhibit stable spurious steady states, as well as spurious stable 
fixed points of period higher than one. In the majority of cases, these occur for values of r 
above the linearized stability limit of the true fixed point. But this is not always the case, as 
is demonstrated in the modified Euler scheme applied to the logistic ODE (3.7) and (3.23) for 
0 < b < 0.5, and the R-K 4 applied to the logistic ODE. For these two methods and ODE 
combinations, stable spurious fixed points occur below the linearized stability limit. In some 
of the instances, these spurious fixed points are outside the interval of the stable and unstable 
fixed points of the ODEs. Others not only lie below the linearized stability limit but also in the 
region between the fixed points of the DEs and so could be very easily achieved in practice. 
For example, in Fig. 3.2b, the modified Euler scheme for the logistic ODE, the linearized 
stability limit of period 1 E is r = 2. But depending on the value of r, two stable fixed points 
of period 1 (one is spurious) can exist at the same time for 0 < r < 1.236. For the R-K 
4 method applied to the logistic DE, one can see from Fig. 3.2d that spurious steady states 
which exist for 2.75 < r < 2.785 are below the linearized stability limit of the l E branch. For 
the modified Euler method applied to du/di = au(l - u)(b - u), it is interesting to see the 
changing behavior of stable spurious steady states as the stable fixed point u* = 6 is varied 
between 0 and 0.5. 

A unified analysis for the above for the standard explicit Runge-Kutta methods is reported in 
Griffiths et al. (1992a). Tables 3.1 - 3.4, taken from Griffiths et al. (1992a), show the true and 
spurious asymptotes of selected schemes. Some entries are marked with an asterisk to indicate 
where stable fixed points are known to exist but no closed analytic form has yet been found. 

Historically, Iserles (1988) was the first to show that while LMMs for solving ODEs possess 
only the fixed points of the original DEs, popular Runge-Kutta methods may exhibit spurious 
fixed points. Iserles et al. (1990) and Hairer et al. (1989) classified and gave guidelines 
and theory on the types of Runge-Kutta methods that do not exhibit spurious period one or 
period two fixed points. Humphries (1991) showed that under appropriate assumptions if stable 
spurious fixed points exist as the time-step approaches zero, then they must either approach 
a true fixed point or become unbounded. Hence repeating the integration with a smaller step 
size will ultimately make spurious behavior apparent. However, convergence in practical 
calculations involves a finite time step A< that is not small as the number of integrations n -> oo 
rather than At — » 0, as n — > oo. The work in Yee et al. (1991), Yee & Sweby (1994, 1995a,b), 
and Lafon & Yee (1991, 1992), Sweby et al. (1990, 1995), Sweby & Yee (1991), and Griffiths 
et al (1992) attempted to provide some of the global asymptotic behavior of time discretizations 
when finite fixed but not extremely small At is used. As mentioned in Section II, constant 
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step size algorithms are relevant to many applications in CFD as well as in other science and 
engineering disciplines. Thus a clear understanding of the dynamics of constant step sizes is 
necessary to isolate some of the sources and cures of numerical uncertainties in computational 
sciences. 


3.3. Bifurcation Diagrams 

This section discusses another method for obtaining the stable fixed point diagrams or 
bifurcation diagrams before illustrating the symbiotic relationship between permissibility of 
spurious steady states and initial data in fixed time step computations. 

“Full” Bifurcation Diagram (“Complete Fixed Point Diagram”): If one obtains the full 
spectrum of these fixed points of any order as a function of the step size, the fixed point 
diagram is sometimes referred to as the “full” bifurcation diagram. In other words, the “full” 
bifurcation diagram exhibits the complete asymptotic solutions of the discretized counterparts 
as a function of the discretized parameter r. In computing the “full” bifurcation diagram, 
searching for the roots and testing for stability of highly complicated nonlinear algebraic 
equations (for fixed points of higher period and/or complex nonlinear DEs combination) can 
be expensive and might lead to inaccuracy. In certain instances, one might be able to obtain 
the bifurcation diagram by some type of continuation method. The most popular one is called 
the pseudo arclength continuation method and was devised by Keller (1977). However, the 
majority of the continuation type methods require known start up solutions for each of the 
main bifurcation branches before one can continue the solution along a specific main branch. 
For problems with complicated bifurcation patterns, the arclength continuation method cannot 
provide the complete bifurcation diagram without the known start up solutions. In fact, it is 
usually not easy to locate even just one solution on each of these branches, especially if spurious 
asymptotes exist. 

Computed Bifurcation Diagram . A numerical approach in obtaining a “computed” bifurca- 
tion diagram (not necessarily the full bifurcation diagram, as explained later) of the resulting 
discretized counterpart is by iterations of the underlying discrete map. In other words, this type 
of computed bifurcation diagram for the one-dimensional discrete map displays the iterated 
solution u n vs. r after iterating the discrete map for a given number of iterations with a chosen 
initial condition for each of the r parameter values. For the figures shown later, with a given 
interval of r and a chosen initial condition, the r axis is divided into 500 equal spaces. In 
each of the computations, the discrete maps were iterated with 600 preiterations (more or less 
depending on the DE and scheme combinations) and the next 200 iterations were over plotted 
in the same diagram for each of the 500 r values. The preiterations are necessary in order for 
the solutions to settle to their asymptotic value. A high number of iterations are overlaid on the 
same plot in order to detect periodic orbits (in this case periods of up to 200) or inva rian t sets. 
The reader is reminded that with this method of computing the bifurcation diagrams, only the 
stable branches are obtained. The domains of the r and u" axes are chosen to coincide with 
the stable fixed point diagrams shown previously. As explained later, even though all of these 
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discrete maps possess periodic solutions of period n for arbitrarily large n and stable chaotic 
solutions, no attempt was made to compute all of the spurious orbits of any order or chaotic 
solutions. The purpose of the present discussion is to show the spurious behavior and these 
computations suffice to serve the purpose. 

j Examples: Figure 3.4 shows the bifurcation diagram of the Euler scheme applied to the logistic 
DE with an initial condition u° > 0. It is of interest to know that in this case the bifurcation 
diagram looks practically the same for any u° > 0. This is due to the fact that no spurious 
fixed points exist for r < 2 and no spurious asymptotes of low period exist for r < 2.627. One 
quickly observes that using the arclength continuation method for this discrete map is the most 
efficient way to obtain its bifurcation diagram. However, this is not the case for other methods 
to be discussed later. Comparing the bifurcation diagram with Fig. 3.2a, one can see that if 
we had computed all of the fixed points of period up to 200 for Fig. 3.2a, the resulting fixed 
point diagram would look the same as the corresponding bifurcation diagram (assuming 600 
iterations of the logistic map are sufficient to obtain the converged stable asymptotes of period 
up to 200 and the chosen initial data are appropriate to cover the basins of all of the periods in 
question). The numeric labeling of the branches in the bifurcation diagram denote their period, 
with only the essential ones labeled for identification purposes. 

The noise appearing on the 1® branch near the bifurcation point r = 2 of the linearized 
stability limit of the fixed point u* = 1 indicates that 600 iterations of the logistic map are not 
sufficient to obtained the converged stable asymptotes. This phenomenon is common to other 
bifurcation points of higher periods as well as the rest of the corresponding bifurcation for the 
studied schemes. See Yee et al. (1991), and Yee & Sweby (1994, 1995a, b) for additional 
details. In fact, the slow convergence of using a time step that is near the linearized stability of 
the scheme (bifurcation point) might be due to this fact. 

It is remarked that the explicit Euler applied to the logistic ODE resulted in the famous 
logistic map. Unlike the underlying logistic ODE, it is well known that the logistic map 
possesses very rich dynamical behavior such as period-doubling (of period 2 n for any positive 
n) cascades resulting in chaos (Feigenbaum, 1978). One can find Fig. 3.4 appearing in 
most of the elementary dynamical systems text books. The exact -values of r for all of the 
period-doubling bifurcation points and chaotic windows (intervals of r) are known and were 
discovered by Feigenbaum in the late 1970’s. Interested readers should consult these elementary 
text books for details. In other words, one can obtain the analytical (exact) behavior of the 
spurious asymptotes and numerical (spurious) chaos of the logistic map. The next section 
explains why using a single initial datum in computing the bifurcation diagrams for schemes 
that exhibit spurious asymptotes does not necessarily coincide with the fixed point diagram (or 
full bifurcation diagram). It is interesting to note that the corresponding bifurcation diagrams of 
the respective discrete maps produced by the remaining studied schemes consist of unions of 
“logistic-map-like” bifurcations and/or “inverted logistic-map- like” bifurcations with similar 
yet slightly complicated period-doubling cascades resulting in chaos. See Section 3.4 for 
additional discussions. 
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Types of Bifurcations . In all of the fixed point diagrams shown previously, the majority of 
the bifurcation phenomena can be divided into three kinds; these are flip, supercritical and 
transcritical bifurcations (Seydel 1988). Figure 3.5 shows examples of these three types of 
bifurcation for the logistic ODE using the modified Euler, improved Euler and R-K 4. Figure 
3.5 also shows a comparison of the stable and unstable fixed points of periods 1 and 2. Although 
the modified Euler and the R-K 4 experience a transcritical bifurcation, they have different 
characteristics. 

For the bifurcation of the first kind, the paths (spurious or otherwise) resemble period 
doubling bifurcations (flip bifurcation) similar to the logistic map. See Fig. 3.2a (forr = 2) and 
Fig. 3.2e (for r = 2) for examples. The second kind is the steady or supercritical bifurcation. It 
occurs, most often, at the main branch 1e, with the spurious paths branching from the correct 
fixed point as it reaches the linearized stability limit, and quite often even bifurcating more 
than once, as r increases still further before the onset of period doubling bifurcations. See 
Fig.3.2c (for r = 2) and Fig. 3.2f (for r = 2) for examples. Using the P-C 3 method to solve 
(3.7) more than one consecutive steady bifurcation occurs before period doubling bifurcations. 
Follow the Is labels on Fig. 3.2f. Although figures are not shown for ODE (3.23) with b = 0.5, 
the improved Euler experiences two consecutive steady bifurcations before period doubling 
bifurcation occurs (see original paper for details). Using the P-C 3 method to solve (3.23), four 
consecutive steady bifurcations occur before period doubling bifurcations. The modified Euler 
and R-K 4 methods, however, experience only one steady bifurcation before period doubling 
bifurcations occur. 

The third kind of bifurcation again occurs most often at the main branch 1 e - The spurious 
paths near the linearized stability limit of 1® experience a transcritical bifurcation. See Fig. 
3.2b (forr = 2), Fig. 3. 2d (forr near 2.75), Fig. 3.2f (forr near 3.4) and Fig. 3.3 (follow the 1 e 
branch) for examples. Notice that the occurrence of transcritical and supercritical bifurcations 
is not limited to the main branch 1 e - See Fig. 3.3 for examples. At the stability limit of the true 
fixed point, only the modified Euler and R-K 4 undergo transcritical bifurcation. 

As can be seen, the occurrence of flip and supercritical bifurcations is more common. In fact, 
most of the bifurcation points shown in previous figures are of these types. The other commonly 
occurring bifurcation phenomenon is the subcritical bifurcation which was not observed in our 
two chosen S(u) functions. With a slight change in the form of our cubic function S(u), a 
subcritical bifurcation can be achieved (Seydel 1988). See elementary text books on bifurcations 
of discrete maps (Seydel 1988) for a discussion of these four types of bifurcation phenomena. 
The consequence of the latter three bifurcation behaviors is that bifurcation diagrams computed 
from a single initial condition u° will appear to have missing sections of spurious branches, 
or even seem to jump between branches. This is entirely due to the existence of spurious 
asymptotes of some period or more than one period, and its dependence on the initial data. 

Slow Convergence Near Bifurcation Points: As discussed previously, the number of itera- 
tions for the computed asymptotes that are near or at the bifurcation point can be orders of 
magnitude higher than away from the bifurcation points. In fact, depending on the type of 
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bifurcation and initial data, one might experience slow convergence using a time step that is 
near the linearized stability of the scheme (bifurcation points of the above four types). See Yee 
& Sweby (1994, 1995a,b) for some examples. In the worst case scenario, if the bifurcation is of 
the transcritical or subcritical type and the time step is within that range, the numerical solution 
can get trapped in a spurious steady state or a spurious limit cycle, causing nonconvergence of 
the numerical solution. Due to the fact that if the spurious steady state or spurious limit cycle 
is very close to the true steady state with a time step that is near the bifurcation point, it is not 
easy to distinguish it from a pure slow convergence issue or a spurious issue. 


3.4. Strong Dependence of Solutions on Initial Data ( Numerical Basins of Attraction) 

Computing Bifurcation Diagrams Using A Single Initial Datum: Figures 3.6a - 3.6c show 
the bifurcation diagram by the modified Euler method for the logistic ODE with three different 
starting initial conditions (I.C.). In contrast to the explicit Euler method as shown in Fig. 3.4, 
none of these diagrams look alike. One can see the influence and the strong dependence of 
the asymptotic solutions on the initial data. For certain initial data and At value combinations, 
spurious dynamics can be avoided. Yet for other combinations, one can never get to the correct 
steady state. In other words, it is possible that for the same At but two different initial data or 
vice versa, the scheme can converge to two different distinct numerical solutions of which one 
or neither of them is the true solution of the underlying DE. Thus in a situation where there is 
no prior information about the exact steady-state solution, and where a time-marching approach 
is used to obtain the steady-state numerical solution when initial data are not known, a stable 
spurious steady-state could be computed and mistaken for the correct steady-state solution. 
Figure 3.6d shows the corresponding ‘ ‘full ’ ’ bifurcation diagram, their earlier stages resembling 
the fixed point diagram 3.2b which showed only fixed points up to period 8. See Yee et al. 
(1991) for an example that overplotting a number of initial data, but not the appropriate ones, 
would not be sufficient to cover all of the essential spurious branches. The strong dependence 
of solutions on initial data is evident from the various examples in which this type of behavior 
is present. It is remarked that if one uses the pseudo arclength continuation type method without 
solving for the roots of the spurious fixed point, one only knows one starting solution (i.e., the 
exact steady states of these two ODEs). The continuation method, in this case, only produces 
the branch of the bifurcation diagram originating from the 1 e branch of the curve. 

Computed Full Bifurcation Diagram : In order to compute a ‘ ‘full’ ’ bifurcation diagram using 
this numerical approach, we must overplot all of the individual bifurcation diagrams of existing 
asymptotes of any period and chaotic attractors obtained by using the entire domain of u values 
as starting initial data. Thus, a better method in numerically approximating the full bifurcation 
diagram is dividing the domain of interest of the tt axis into equal increments and using these u 
values as initial data. The “full” bifurcation diagram is obtained by simply overplotting all of 
these individual diagrams on one. Figs. 3.7 and 3.8 show the “full” bifurcation diagrams for the 
corresponding fixed point diagrams shown previously. Note that the full bifurcation computed 
this way might miss some of the windows of bifurcations that occur inside the intervals of the 
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adjacent r and/or the initial data values. 

It is noted that for the cases when one knows the bifurcation pattern of a specific discrete map, 
the actual number of the initial data points used for that full bifurcation diagram computations 
do not have to completely cover the entire domain of u as long as these initial data cover all 
of the basins of attraction of the asymptotes (which initial data lead to which asymptotes). See 
next section for definition and discussion. That is, at least one initial data point is used from 
each of the basins of the asymptotes. No attempt has been made to compute the complete full 
bifurcation diagram, since this is very costly and involves a complete picture of the existing 
asymptotes of any period and chaotic attractors for the domain of interest in question. See 
remarks in Section 3.2 on computed bifurcation diagrams for explanation. Here, we use the 
term “full bifurcation diagram” to mean “computed bifurcation diagram with sufficient initial 
data to cover the essential lower order periods”. Without loss of generality, from here on we 
use the word bifurcation diagram to mean the computed (and approximated) full bifurcation 
diagram. 

From Figs. 3.7 and 3.8, one can conclude that all of the studied explicit methods undergo 
“period doubling bifurcations” leading to the “logistic-map-type bifurcations”. The term 
“logistic-map-type bifurcations” here means that the behavior and shape of the bifurcations 
resemble the logistic map as shown in Fig. 3.4. The ranges of the r values in which logistic- 
map-type bifurcations occur are not restricted to the 1 e branch of the bifurcation diagram. The 
birth of the logistic-map-type bifurcations can occur below or beyond the linearized stability 
limit of the true steady state of the governing equation. To aid the readers. Figs. 3.7 and 3.8 
indicate the major stable fixed points of periods up to 4. Basically (not restricted to), each 
of the Is branches of the bifurcations are logistic-map-type. For example, the modified Euler 
experiences two period doubling bifurcations for the ODE (3.7). For the ODE (3.23), the 
modified Euler experiences at least two to three period doubling bifurcations, depending on the 
b values. For other methods, the situation is slightly more complicated. 

Besides the regular logistic-map-type bifurcations, some of these methods undergo the so 
called “inverted logistic-map- type” of period doubling bifurcations. The shape of this type 
of bifurcation resembles the reverse image of Fig. 3.4. See, for examples. Fig. 3.7b for 
1.62 < r < 1.67 and Fig. 3.7c for 3.15 < r < 3 . 3 . 

The above example explains the role of initial data in the generation of spurious steady- 
state numerical solutions, stable and unstable spurious numerical chaos and other asymptotes. 
Section 3.5 illustrates the role of initial data in the permissibility of stabilizing unstable 
steady states of the governing equation and the introduction of stable and unstable spurious 
numerical chaos and other asymptotes by implicit LMMs. Next, how basins of attraction can 
complement the bifurcation diagrams in gaining more detailed global asymptotic behavior of 
time discretizations for nonlinear DEs will be illustrated. 

“Exact” Basin of Attraction vs. “Numerical” Basin of Attraction: The basin of attraction 
of an asymptote (for the DEs or their discretized counterparts) is a set of all initial data 
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asymptotically approaching that asymptote. We use the terms “exact” basin of attraction and 
‘ ‘numerical” basin of attraction to mean the basin of attraction of the DE and basin of attraction 
of the underlying discretized counterpart. Although all of the numerical basins of attraction 
shown later are obtained numerically, the term ‘ ‘numerical basin of attraction” here pertains to 
the true basins of attraction of all of the asymptotes for the underlying discrete map. 

For the logistic ODE, the “exact” basin of attraction for the only stable fixed point u* = 1 
is the entire positive plane for all values of a > 0. The basin of attraction for the ODE (3.23) for 
the stable fixed point u* = b is 0 < u < 1 for all a > 0, 0 < b < 1. However, the situation for 
the corresponding numerical basins of attraction for the various schemes is more complicated 
since each asymptote, stable or unstable, spurious or otherwise, possesses its own numerical 
basin of attraction. Intuitively, in the presence of stable and unstable spurious asymptotes, 
the basins of the true stable and unstable steady states (and asymptotes if they exist) can be 
separated by the numerical basins of attraction of the stable and unstable spurious asymptotes. 
Consequently, what is part of the true basin of a particular fixed point of the governing equation 
might become part of the basin of the spurious asymptotes. For implicit methods that can 
stabilize unstable steady states of the governing equation (to be discussed later), the basin 
of attraction associated with the particular stabilizing steady state can consist of the union of 
part of basins from other true asymptotes. In other words, the allowable initial data of the 
numerical method associated with a particular asymptote of the DE can experience enlargement 
or contraction, and can become null or consist of a union of disjoint regions. These regions can 
be fractal like. Therefore, keeping track of which initial data lead to which asymptotes (exact or 
spurious) of the underlying discrete maps becomes more complicated as the number of spurious 
asymptotes increases. 

On the other hand, the computed bifurcation diagram cannot distinguish the types of 
bifurcation and the periodicity of the spurious fixed points of any order. With the numerical 
basins of attraction and their respective bifurcation diagrams superimposed on the same plot, 
the type of bifurcation and which initial data lead to which stable asymptotes become apparent. 
In order to obtain the corresponding numerical basins of attraction for the schemes discussed 
above, one immediately realizes that, in most cases, a numerical approach is the only recourse 
until more theoretical tools for searching for the basin boundaries of general discrete maps 
become available. We would like to add that there are isolated theories or approximate methods 
to locate some basin boundaries for simple discrete maps or special classes of discrete maps. 
Even in this case, these methods are neither practical nor are there fixed guidelines for the actual 
implementation of discrete maps for more complex ones of similar type. See Hsu (1987) for an 
approximate method and Friedman ( 1 995) for numerical algorithms which compute connecting 
orbits. 

Computing Num erical Batins of Attraction on the Connection Machine . The nature of 
this type of computation, especially for systems of nonlinear ODEs or PDEs, requires the 
performance of a very large number of simulations with different initial data; this can be 
achieved efficiently by the use of the highly parallel Connection Machine (CM-2 or CM- 5) 
or the IBM SP2 machine whereby each processor could represent a single initial datum and 
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thereby all the computations can be done in parallel to produce detailed global stability behavior 
and the resulting basins of attraction. With the aid of highly parallel Connection Machines, we 
were able to detect a wealth of the detailed nonlinear behavior of these schemes for systems of 
ODEs and PDEs which would have been overlooked had isolated initial data been chosen on 
the Cray or other serial or vector machine. 

Figure (3.9) shows the bifurcation diagram with numerical basins of attraction superimposed 
for the logistic ODE for four Runge-Kutta methods. (Even though one might not need to use a 
highly parallel machine to compute the basins of attraction for the scalar ODE, nevertheless this 
figure was computed on the CM-5 requiring only a few minutes for each plot). The major work 
on the CM-5 coding is on the efficient handling of data for plotting. The same plots would have 
required at least an order of magnitude more CPU time on a serial or vector machine. To obtain 
a bifurcation diagram with numerical basins of attraction superimposed using the CM-5, the 
preselected domain of initial data and the preselected range of the At parameter are divided into 
512 equal increments. For the bifurcation part of the computations, the discretized equations 
are, in general, preiterated 3000 steps for each initial datum and At before the next 1000 
iterations (more or less depending on the problem and scheme) are plotted. The preiterations 
are necessary in order for the trajectories to settle to their asymptotic value. The high number 
of iterations are plotted (overlaid on the same plot) to detect periodic solutions. The bifurcation 
curves appear on the figures as white curves, white dots and dense white dots. While computing 
the bifurcation diagrams it is possible to overlay basins of attraction for each value of A t used. 
For the numerical basins of attraction part of the computations, with each value of At used, 
we keep track where each initial datum asymptotically approaches, and color code each basin 
(appearing as a multi-color vertical line) according to the individual asymptotes. Black regions 
denote basins of divergent solutions. While efforts were made to match color coding associated 
with a particular asymptote of adjacent multi-color vertical lines on the bifurcation diagram 
(i.e., from one At to the next), it was not always practical or possible. Care must therefore be 
taken when interpreting these overlays. In an idealized situation it is best if we also know the 
critical value of At for the onset of unstable spurious asymptotes. However, with the current 
method of detecting the bifurcation curve only the stable ones are detected. For example, a 
steady bifurcation would break the domain immediately after the bifurcation point into two 
different color domains immediately after the bifurcation point. 

Examples : Any initial data residing in the green region in Fig. 3.9 for the modified Euler 
method belongs to the numerical basin of attraction of the spurious (stable) branch emanating 
from u = 3 and r = 1 . Thus, if the initial data is inside the green region, the solution can never 
converge to the exact steady state using even a small fixed but finite At (all below the linearized 
stability limit of the scheme). Although not shown, we have computed the bifurcation diagram 
with wider ranges of r and initial data. In fact, the green region actually extends upward as 
r decreases below 1. Thus for a small range of r values very near zero, the entire domain is 
divided into two basins (same as the ODE). As r increases, the domain can divide into more 
than two basins (instead of two for the ODE). But of course higher period spurious fixed points 
exist for other ranges of r and more basins are created within the same u domain. For r near 


28 



2 (i.e., near the linearized stability limit of the true steady state u* = 1, the bifurcation is 
transcritical. Using an r value slightly bigger than 2 will lead to the spurious steady state until 
r increases beyond 3.236. Consequently, there are large ranges of r below and beyond the 
linearized stability limit of the true steady state for which spurious dynamics occur. Observe 
the size and shape of these basins as r varies. The majority of these basin boundaries are fractal 
like. 

A similar situation exists for the R-K 4 method (Fig. 3.9), except that now the numerical 
basins of attraction of the spurious fixed points occur very near the linearized stability limit 
of the scheme, with a small portion occurring below the linearized stability limit. Although 
both the modified Euler and the R-K 4 methods experience a transcritical bifurcation for the 
logistic ODE. The transcritical bifurcation for the R-K 4 is more interesting. See Fig. 3.5 for the 
distinction between the two transcritical bifurcations. In contrast to the improved Euler method, 
the green region represents the numerical basin of attraction of the upper spurious branch of 
the transcritical bifurcation. The bifurcation curve directly below it with the corresponding red 
portion is the basin of the other spurious branch. With this way of color coding the basins of 
attraction, one can readily see (from the plots) that for the logistic ODE (3.7), the improved 
Euler method experience one steady bifurcations before period doubling bifurcation occurs. 
The modified Euler and R-K 4 methods, however, experience a transcritical bifurcation before 
period doubling bifurcations occur. The Kutta method experiences period doubling bifurcation 
at the linearized stability limit. One way to detect the steady bifurcation from these plots is to 
look for a separate color associated with each branch of the associated bifurcation. A similar 
interpretation holds for certain type of transcritical bifurcation. See R-K 4 in Fig. 3.9. One 
way to detect the period doubling bifurcation from these plots is to look for no change in 
colors associated with each branch of the associated bifurcation. For subcritical bifurcation, it 
is slightly more complicated. 

The above discussion shows the interplay between initial data, step size and permissibility of 
spurious asymptotes. It indicates that it is not just the occurrence of stable spurious numerical 
solutions that causes difficulty. Indeed such cases may be easier to detect. These spurious 
features of the discretizations often occur but are unstable; i.e., they do not appear as an actual 
(spurious) solution because one usually cannot obtain an unstable asymptotic solution by mere 
forward time integration. However, far from being benign, they can have severe detrimental 
effects on the allowable initial data of the true solution for the particular method, hence causing 
slow convergence or possibly even nonconvergence from a given set of initial data, even though 
the data might be physically relevant. 

Due to space limitations, interested readers are referred to Yee & Sweby (1994, 1995a) 
for results for four 2x2 systems of nonlinear model ODEs. Classification of fixed points 
of systems of equations are more involved than the scalar case. See elementary dynamical 
systems text books for details. The corresponding bifurcation diagrams and numerical basins 
of attraction of these schemes are even more involved. New phenomena exist that are absent 
from the scalar case. For example, in the presence of multiple steady states, even for explicit 
methods, depending on the step size and initial data combination, the associated numerical 
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basins of attraction for a true steady state might experience an enlargement of their basins at the 
expense of a contraction of the other asymptotes. For the scalar case, this is only possible for 
superstable implicit methods (see next section). See Yee & Sweby (1994, 1995a,b) for more 
details. Another example is that the fixed point can change type as well as stability (e.g., from 
a saddle point to a stable or unstable node). It is remarked that for systems beyond 3 x 3, it is 
impossible to conduct a detailed analysis as shown above. 

3.5. Global Asymptotic Behavior of Superstable Implicit LMMs 

This section reviews the superstable property of some implicit LMMs and summarizes their 
global asymptotic nonlinear behavior using the dynamical approach. Recall that the underlying 
discrete maps from using LMMs are linear in the discretized parameter A t and they are exempt 
from spurious steady states. As can be seen later, the combination of implicit LMMs and the 
superstable property produce asymptotic behavior that is very different from schemes that were 
studied in the previous section. One distinct property of these types of schemes is that they 
can stabilize unstable steady states of the governing equation. Another property is that the 
numerical basins of attraction of the stable steady state can include regions of the basins of the 
unstable steady states of the governing equation. 


3.5.1. Super-stability Property 

Dahlquist et al. (1982) first defined super-stability in ODE solvers to mean the region 
of numerical stability that encloses regions of physical instability of the true solution of the 
ODE. Dieci & Estep (1991) subsequently gave a broader definition as one in which an ODE 
solver does not detect that the underlying solution is physically unstable. They observed that 
super-stability can occur also when the ODE solver is not super-stable in terms of Dahlquist 
et al. They concluded that the key factor which determines the occurrence of super-stability 
is the iterative solution process for the nonlinear algebraic equations. They indicated that the 
iterative solution process has its own dynamics, which might be in conflict (as for Newton’s 
method) with the dynamics of the problem. They also indicated that super-stability can arise 
because of this fact. Their viewpoint is that the numerical scheme and the methods for solving 
the nonlinear algebraic equations should be considered as a unit. Neglecting this latter aspect 
and basing step size selection purely on accuracy considerations leads to faulty analysis. They 
believe that error control strategies for stiff initial value problems ought to be redesigned to take 
into account stability information of the continuous problem. 

In Yee & Sweby (1994, 1995a,b) we exploited the global asymptotic behavior of some of 
the superstable implicit LMMs for constant step sizes. We concentrated on four of the typical 
unconditionally stable implicit LMMs. These are backward Euler, trapezoidal rule, midpoint 
implicit, and three-point backward differentiation (BDF), each with noniterative (linearized 
implicit), simple, Newton and modified Newton iterative procedures for solving the resulting 
nonlinear algebraic equations. A semi-implicit predictor method of Gresho (1996) also was 
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investigated. 

We believe that some of the phenomena observed in our study were not observed by Died 
& Estep (1991) or by Iserles (1988). Based on our study, we now give a loose definition of an 
implicit time discretization as having a super-stability property if, within the linearized stability 
limit for a combination of initial data and time step (fixed) or starting time step (using standard 
variable time step control based on accuracy requirements), the scheme stabilizes unstable 
steady states of the governing equations in addition to having the property of Died and Estep. 
The definition includes the procedures in solving the nonlinear algebraic equations. This loose 
definition fits the behavior that was observed in Yee & Sweby (1994, 1995a,b), while at the 
same time it fits the framework of dynamics of numerics in time-marching approaches in CFD. 
This is not a re-definition of Died and Estep, but rather a clarification of their definition when 
asymptotic numerical solutions of the governing equations are desired. In this case, superstable 
schemes might have the property of a numerical basin of attraction of the true steady state that 
is larger than the underlying exact basin of attraction. As can be seen in Yee & Sweby (1994, 
1995a,b), the trapezoidal method is more likely to exhibit this property than the other three 
LMMs. The stabilization of unstable steady states by LMMs was also observed by Salas et al. 
(1986), Embid et al. (1984), Burton & Sweby (1995) and Poliashenko (1995). Section 4.2.2 
gives a summary of the work of Burton and Sweby. 

3.5.2. Implicit LMMs 


The four LMMs and a semi-implicit predictor method considered are 

Implicit Euler Method 


u n+1 = tx n + AiS n+1 , 

(3.24) 

Trapezoidal Method 


u n+1 = u n + lr(S n + S n+1 ), 

(3.25) 

Midpoint Implicit Method 


U -‘-u" + r5( ( “" +, 2 +Un) ), 

(3.26) 

3-Level Backward Differentiation Formula (BDF) 


u n+1 = u" + ^rS n+1 + ^(u n - ti" -1 ), 

O O 

(3.27) 


Semi-Implicit Predictor Method of Gresho (Gresho 1 996). 

The semi-implicit predictor method is the same as (3.24) but with an added predictor step using 
the explicit Euler before the implicit step (3.24) to make the final scheme second-order accurate. 
The four methods of solving the resulting nonlinear algebraic equations are as follows. 
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Linearization (a noniterative procedure) is achieved by expanding S n+1 as 5 n + 
J r (u n )(u n+1 — u n ), where J(u n ) is the Jacobian of S. 

Simple Iteration is the process where, given a scheme of the form u n+1 = G(u n , u n+1 ), 
we perform the iteration 


= (3-28) 

where «q + 1 = u n and “(i/)” indicates the iteration index. The iteration is continued either 
until some tolerance between iterates is achieved or a limiting number of iterations has been 
performed. In all of our computations the tolerance “tol” is set as Hu"* 1 - u”J t ^ , 1) || < tol and 
the maximum number of iterations is 15. The major drawback with simple iteration is that for 
guaranteed convergence the iteration must be a contraction; i.e. 


||G(u V) - G(«”, tr)|| < a||t, - w ||, (3.29) 

where a < 1. Whether or not the iteration is a contraction at the fixed points will influence 
the stability of that fixed point, overriding the stability of the implicit scheme. Away from the 
fixed points the influence will be on the basins of attraction. In other words, “implicit method 
+ simple iteration” behaves like explicit methods. As can be seen in Yee & Sweby (1994, 
1995a,b), numerical results illustrate this limitation in terms of basins of attraction as well. One 
advantage of the “implicit method + simple iteration” over non-LMM explicit methods is that 
spurious steady states cannot occur. Due to this fact, results using simple iteration will not be 
presented here. 

Newton Iteration for the implicit schemes is of the form 


u 


n+l _ n+1 

(^+i) - V) 




(3.30) 


where «o +1 = u ”- The differentiation is with respect to the second argument and the scheme 
has been written in the form (for two- level schemes) F(u n ,u n+1 ) = 0. 


Modified Newton iteration is the same as (3.30) except it uses a frozen Jacobian F'(u n , u"). 
The same tolerance and maximum number of iterations used for the simple iteration are also 
used for the Newton and modified Newton iterations. In all of the computations, the starting 
scheme for the 3-level BDF is the linearized implicit Euler. 


We also considered two variable time step control methods. The first one is “implicit Euler 
+ Newton iteration with local truncation error control” 


•S* = »’ 


u n+l _ n+1 

V+i) ~ V) + 


J-Arj(u^) 


tift 1 - - A t n S{u”+ l ) 


1 ,... 

(3.31a) 
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with 


A t n = O.gAr-^toU/Hu" - a”" 1 - A< n - 1 5(« n )||} 1/3 , (3.31b) 

where the (n + l)th step is rejected if ||u n — u n_1 — A< n_1 5’(«")|| > 2toli. In this case, we 
set A t n = A < n_1 . The value “toll” is a prescribed tolerance and the norm is an infinity norm. 
The second one is the popular “ode23” method 

fc, = S(u n ) 
k 2 = S(u n + At n k t ) 
k s = S(u n + At n {kt + k 2 )/4) 
u n+1 = u” + A< n (*i + k 2 + 4*,)/6 
Au n+1 = At n (ki + k 2 — 2ki)/Z (3.32a) 

with 

A ‘" = 0 - 9A ‘""‘\/5S- (3 - 32b) 

where the (n + l)th step is rejected if ||Au T,+1 || > toll max{l, ||u n+1 ||}. In that case, we set 
A/ n_1 = A < n . Again, “toll” is a prescribed tolerance and the norms are infinity norms. We 
also use the “straight Newton’’ method to obtain the numerical solutions of 5(u) = 0, which 
is the one-step Newton iteration of the implicit Euler method of (3.30). 

We mapped out the bifurcation diagrams and numerical basins of attraction for these 
five schemes as a function of the time step for different nonlinear model equations with 
known analytical solutions (scalar and 2x2 systems of autonomous nonlinear ODEs). The 
computations were performed on the CM-5. In general, we preiterated the discretized equations 
3000 - 5000 steps before the next 3000 iterations are plotted. Comparing the results with the 
explicit methods, it was found that aside from exhibiting spurious asymptotes, all of the four 
implicit LMMs can change the type and stability (unstable to stable) of the steady states of 
the differential equations (DEs). See the next section for some examples. They also exhibit a 
drastic distortion but less shrinkage of the basin of attraction of the true solution than standard 
non-LMM explicit methods. In some cases with smaller At, the implicit LMMs exhibit 
enlargement of the basins of attraction of the true solution. Overall, the numerical basins of 
attraction of the linearized implicit procedure mimics more closely the basins of attraction of 
the continuum than the studied iterative implicit procedures for the four implicit LMMs. In 
general the numerical basins of attraction bear no resemblance to the exact basins of attraction. 
The size can increase or decrease depending on the time step. Also, the possible existence of the 
largest numerical basin of attraction that is larger than the exact one does not necessarily occur 
when the time step is the smallest. Although unconditionally stable implicit methods allow a 
theoretically large time step At, the numerical basins of attraction for large At sometimes are 
so fragmented and/or so small that the safe (or practical) choice of At is only slightly larger or 
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comparable to the stability limit of standard explicit methods (but with larger numerical basins 
of attraction than the explicit method counterparts). In general, if one uses a At that is a fraction 
of the stability limit, one has a higher chance of convergence to the correct asymptote than when 
one uses standard explicit methods. Comparing the results of Yee & Sweby (1994, 1995b) 
with Yee & Sweby (1995a), the implication is that unconditionally stable implicit methods are, 
in general, safer to use and have larger numerical basins of attraction than explicit methods. 
However, one cannot use too large a time step since the numerical basins of attraction can be 
so small and/or fragmented that the initial data has to be very close to the exact solution for 
convergence. This knowledge improved the understanding of the basic ingredients needed for 
a time-marching method using constant step sizes to have a rapid and guaranteed convergence 
to the correct steady states. 

Numerical experiments performed on the two variable time step control methods also 
indicated that, although variable time step controls are not foolproof, they might alleviate 
the spurious dynamics most of the time. One shortcoming is that in order to avoid spurious 
dynamics, the required size of At is impractical to use in CFD, especially for the explicit 
ode23 method. The next section shows some representative global asymptotic behavior of these 
implicit LMMs. 


3.5.3. Numerical Examples 

Selected results in the form of bifurcation diagrams and basins of attraction are shown in 
Figs. 3.10-3.16 for the logistic ODE model (3.7) and Figs. 3.17-3.19 for the second model 
(3.23). In Figs. 3.10, 3.12, 3.14, 3.16-3.19, the left diagrams show the bifurcation diagrams and 
the right diagrams show the bifurcation diagrams with basins of attraction superimposed on the 
same plot. However, Figs. 3.1 1, 3.13 and 3.15 show only the bifurcation diagrams with basins 
of attraction superimposed on the same plot. In all of Figs. 3.10-3.19, the abscissa is r = a At. 

Note that the preselected regions of At and the initial data were determined by examining 
a wide range of At and initial data. In most cases, we examined At from close to zero up to 
one million. What is shown in these figures represents some of the At and initial data ranges 
that are most interesting. Due to this fact, the Ai and initial data ranges shown are different 
from one method to another for the same model problem. The streaks on some of the plots are 
either due to the non-settling of the solutions within the prescribed number of iterations or the 
existence of small isolated regions of spurious asymptotes. Due to the high cost of computation, 
no further attempts were made to refine their detailed behavior since our purpose was to show 
how, in general, the different numerical methods behave in the context of nonlinear dynamics. 

Due to the method of tracking the basins of attraction the color coding of the basin of 
attraction associated with a particular asymptote might vary from one vertical line to the next 
vertical line (i.e, from one At to the next). This is the case for Figs. 3.10, 3.1 1, 3.14, 3.18 and 
3.19. For example in Fig. 3.10 the basin of attraction (as a function of At) for the steady state 
u = 1 is the red region before the appearance of the light blue strip. After the appearance of the 
blue strip, it is the region above the envelope that separates the green and red regions. When in 
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doubt, use the bifurcation diagram as a guide and identify the r value where the sudden birth of 
stable spurious asymptotes occured. Incidently, for this particular discrete map, this envelope 
and the critical value r which undergoes stabilization of the fixed point u = 0 can be obtained 
analytically. Independently, Arriola (1993) derived the analytical form of the envelope. 

The midpoint implicit method behaves similarly to the trapezoidal method. In fact their 
linearized forms are identical. From here on, the midpoint implicit method is not discussed. 

As mentioned before, for unconditionally stable LMMs, the scheme should not experience 
any steady bifurcation from the stable branch because unconditionally stable LMMs preserve 
the stability of the stable steady states. This rule does not apply to unstable steady states 
using super-stable LMMs. Before stabilization of the unstable steady state, super-stable LMMs 
typically undergo “inverted period doubling bifurcations” or the “inverted logistic-map- type 
bifurcations” (or crisis in terms of Grebogi et al. 1983). See Fig. 3.10 for 0.2 < u n < —1 for 
an example. For the ODE (3.23), all of the implicit methods experience at least two inverted 
logistic-map-type bifurcations. See Figs. 3.17-3.19. From these figures, we can observe that all 
of the studied implicit methods can introduce stable spurious chaos since a logistic-map-type of 
spurious bifurcations occur. 

Figures 3.10, 3.11, 3.14, 3.16-3.19 show other situations where the rest of three implicit 
LMMs and the Gresho method stabilize the unstable steady states of the ODEs. It appears that 
the onset of stabilization of the unstable steady states arises in two ways. One way is the birth 
of stable spurious asymptotes or stable spurious chaos in the form of an inverted logistic map. 
The second way is the birth of unstable spurious asymptotes (fixed points other than period 
one) leading to the onset of stabilization of unstable steady states. Although the two ways of 
stabilization of the unstable steady states are similar, their corresponding basins of attraction 
are very different. See Figs. 3.10, 3.14, 3.16-3.19. 

The critical value of A t c for the onset of stabilization is not very large. It is problem 
dependent, method dependent and also procedure dependent (the various ways of solving the 
nonlinear algebraic equations). In most cases, the value of the A t c is comparable to or smaller 
than the equivalent of the stability limit of standard explicit Runge-Kutta methods. It is not 
uncommon for the underlying basins of attraction to be larger than the exact basin of attraction 
for At < A t e . 

Among the three procedures, the linearized noniterative forms have a higher tendency to 
stabilize unstable steady states. See Figs. 3.10, 3.14, 3.16-3.19. Here the word “procedure” 
excludes the simple iteration method. Among the three LMMs, the trapezoidal method is the 
least likely to stabilize unstable steady states, but the corresponding basins of attraction can be 
very small and more fragmented than for the other two LMMs. Also, the A t c for stabilization is 
bigger than for the other two LMMs. For a particular LMM not all three procedures for solving 
the nonlinear algebraic equation necessarily stabilize the unstable steady states (see Figs. 3.14 
and 3.15). But, if they do, the pattern or the method for the onset of the stabilization does not 
have to be the same (see Figs. 3.10 and 3.1 1) but the value of A t c is the same for all models 
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and methods studied. 


For the case of the semi-implicit predictor method, the onset of the stabilization can be 
accompanied by the birth of other spurious asymptotes (other than steady state). See Fig. 3.16 
for r > 2. It is fascinating to see how complicated the basins of attraction are which compose 
the many disjoint and fractal like regions. Similar behavior is also observed for the “implicit 
Euler + modified Newton” but is less pronounced than the semi-implicit predictor method. See 
Figs. 3.1 1 and 3.16. 

Compared with the three implicit LMMs, and independent of the method of solving the 
nonlinear algebraic equation, the Gresho method exhibits the smallest basin of attraction 
(compared with the exact basin of attraction of the stable steady state) and is more fragmented 
for A t < A t e . Aside from the efficiency issue, the implication is that a higher order accuracy 
scheme might not be as desirable for the time-marching approach. 

Since straight Newton is just a one step “implicit Euler + Newton”, its basins of attraction 
(for A t larger than explicit Runge-Kutta) are almost the same even with more than one step 
iterations. Studies indicated that contrary to popular belief the initial data using the straight 
Newton method may not have to be close to the exact solution for convergence. Straight Newton 
exhibits stable and unstable spurious asymptotes. Initial data can be reasonably removed from 
the asymptotic values and still be in the basin of attraction. However, the basins can be 
fragmented even though the corresponding exact basins of attraction are single closed domains. 
See Fig. 6.25 of Yee & S weby ( 1 995a). The cause of nonconvergence may just as readily be due 
to the fact that the numerical basins of attraction are fragmented. If one uses a time step slightly 
bigger than the stability limit of standard explicit methods for the three LMMs, straight Newton 
can have similar or better performance. In fact, using a large A< with the linearized implicit 
Euler method or the implicit Euler + Newton procedure has the same chance of obtaining the 
correct steady state as the straight Newton method if the initial data are not known or arbitrary 
initial data are taken. 

A consequence of all of the observed behavior is that part or all of the flow pattern can 
change topology as the discretized parameter is varied. An implication is that the numerics 
might predict, for example, a nonphysical reattachment flow. Thus even though LMMs preserve 
the same number (but not the same types or stability) of fixed points as the underlying DEs, the 
numerical basins of attraction of LMMs do not coincide with the exact basins of attraction of 
the DEs even for small At. Some of the dynamics of the LMMs observed in our study can be 
used to explain the root of why one cannot achieve the theoretical linearized stability limit of 
the typical implicit LMMs in practice when solving strongly nonlinear DEs (e.g., in CFD). 

3.6. Does Error Control Suppress Spuriosity? 

The previous sections discussed mainly the spurious behavior of long time integrations of 
initial value problems of nonlinear ODE solvers for constant step sizes. The use of adaptive 
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step size based on local error control for implicit methods was studied by Dieci & Estep 
(1991). Dieci and Estep concluded that for superstable LMMs with local step size error control 
and depending on the method of solving the resulting nonlinear algebraic equations, spurious 
behavior can occur. Our preliminary study on the two variable step size control methods 
discussed in the previous section indicated that one shortcoming is that the size of At needed 
to avoid spurious dynamics is impractical (too small) to use, especially for the ode23 method. 
Theoretical studies on the adaptive explicit Runge-Kutta method for long time integration have 
been gaining more attention recently. Recent work by Stuart (1994, 1995), Humphries (1992), 
Higham and Stuart (1995) and Aves et al. (1995) showed that local error control offers benefits 
for long-term computations with certain problems and methods. Aves et al. addressed the heart 
of the question of whether local error control confers global properties of steady states of the 
IVP of autonomous ODEs using adaptive Runge-Kutta type methods. 

Aves et al.’s work is concerned with long term behavior and global quantities of general 
explicit Runge-Kutta methods with step size control for autonomous ODEs. They believed that 
the limit t n — > oo is more relevant than the limit of the variable step sizes h n —* 0. They 
studied spurious fixed points that persist for arbitrarily small error tolerances r. This type 
of adaptive Runge-Kutta method usually consists of a primary and secondary Runge-Kutta 
methods of different order. Their main result is positive. When standard local error control is 
used, the chance of encountering spuriosity is extremely small. For general systems of ODEs, 
the constraints imposed by the error control criterion make spuriosity extremely unlikely. 
For scalar problems, however, the mechanism by which the algorithm succeeds is indirect — 
spurious fixed points are not removed, but those that exist are forced by the step size selection 
mechanism to be locally repelling (with the relevant eigenvalues behaving like 0(l/r)). 

To be more precise, adaptive time stepping with Runge-Kutta methods involves a pair of 
Runge-Kutta formulae and a tolerance parameter "r", which is usually small. See for example 
the “ode23” method (3.32). Hence a spurious fixed point of the adaptive procedure requires: 

1) A spurious fixed point common to both methods must exist. This is usually easy to achieve 
since the bifurcation diagrams of individual Runge-Kutta methods have so many branches. 

2) This spurious fixed point must be stable. This is much more difficult to achieve - essentially 
since the bifurcation curves for the two methods must intersect tangentially; otherwise there 
will be an eigenvalue of the Jacobian of 0(l/r). 

They then show that the probability of 2) occurring is zero. However, for a given pair one 
can generally construct functions for which it holds (generally stability will only hold for r > 
lower bound). In any event, the basin of attraction of this spurious fixed point will only be 
0(t). These results were derived for scalar problems. 

In the worst scenario problem classes exist where, for arbitrary r, stable spurious fixed points 
persist with significant basins of attraction. They derived a technique for constructing ODEs for 
which an adaptive explicit Runge-Kutta method will behave badly. They showed that this can 
be accomplished using a locally piecewise constant function 5(u). Since the disjoint pieces can 
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be connected in any manner, S can be made arbitrarily smooth. Hence, smoothness of S alone is 
not sufficient to guarantee that spurious behavior will be eliminated. These examples highlight 
the worst-case behavior of adaptive explicit Runge-Kutta methods. They also mentioned that 
they can construct similar examples involving systems. However, these types of examples are 
somewhat contrived. 

Griffiths is currently working on the application to hyperbolic PDEs. Preliminary results 
(private communication with David Griffiths ( 1 996)) showed that it is by no means clear at the 
moment whether stable spurious solutions may be eliminated. The difference is that, unlike 
physical problems governed by nonlinear ODEs, nonlinear PDEs may have wave-like solutions 
rather than fixed points due to the spatial derivatives. 


3.7. Dynamics of Numerics of a Reaction-Convection Model 

This section further studies the dynamics of selected finite difference methods in the 
framework of a scalar model reaction-convection PDE (LeVeque & Yee 1990) and investigates 
the possible connection of incorrect propagation speeds of discontinuities with the existence 
of some stable spurious steady- state numerical solutions. The effect of spatial as well as time 
discretizations on the existence and stability of spurious steady-state solutions will be discussed. 
This is a summary of the work of Lafon & Yee (1991, 1992). 

A nonlinear reaction-convection model equation in which the exact solution of the governing 
equations are known (LeVeque & Yee 1990) is considered. The model considered in LeVeque 
and Yee is 


u t + au. = aS(u) 
u(z,0) = u°(z) 

where a and a are parameters, and S(u) = — u(l 
by 


0 < x < 1 (3.33a) 

(3.33b) 

u)(2 — u). The boundary condition given 


u(0 ,t) = «o t > 0 
or, periodic boundary condition given by 


(3.33c) 


u(0,<) = u(l,<) <>0 (3.33d) 

is used to complete the system. 

The exact steady-state solutions u* of the continuum PDE (3.33) can be obtained by 
integration by parts of adu* jdx = aS(u *) which yields 


ax 

(- c = 

a 



= log 


K(») - il 1 

VM*)(2-u*(z))|J’ 


(3.34) 
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where c is the integration constant. 


In the case where the boundary condition u(0, t) = u 0 , there is a unique steady state and its 
value is determined by u 0 . If u 0 is a root of 5, (i.e., u 0 = 0, 1, or 2), then the exact steady 
state is constant in x and is equal to u 0 . But if u 0 is not a root of S, then the exact steady state 
satisfies 


x = 


-log 

a 


u*(se) - 1 

/ u 0 (u 0 - 2) 

. ("o - 1) ^ 

j u*(*)(u*(*) — 2) 


(3.35) 


Although, the domain is confined to 0 < * < 1, the steady-state solution is defined for 
0 < * < oo. The limit of u*(z) is 0 as x — ► oo for — oo < u 0 < 0 or 0 < u 0 < 1. The limit of 
u*{x) is 2 if 1 < u 0 < 2, or 2 < u 0 < oo. One can show that this exact steady-state solution is 
stable. 


In the case of the periodic boundary condition where u(0,<) = u(l ,<) and «*(0) is not a root 
of S, it can be shown that there exist three exact steady-state solutions; namely «*(*) = 0, 1, 
or 2. One can also show that u* = 0 and u* = 2 are stable while u* = 1 is unstable. 

Denote the basin of attraction for the steady state u* by BA(u*). Then it is obvious that 
5A(0), the basin of attraction for the steady-state solution u* = 0, is the set of initial data 
u°(x) < 1 for all real values of x. In mathematical notation 

BA(0) = {«° : ti°(») <1 V*}. (3.36) 

Similarly, the basin of attraction for the steady state u* = 2 is 

BA( 2) = {u° : u°(«e) > 1 V*}. (3.37) 

Later we contrast these exact basins of attraction with the numerical basins of attraction SA(0) 
and BA(2) for the various schemes under discussion. 

For the numerical methods, semi-discrete (method of lines) finite difference methods (FDMs) 
and implicit treatment of the source terms (semi-implicit) with noniterative linearization using 
a characteristic form are considered. 


Spatial Discretizations: Let u^(f) represent an approximation to u(*_,-,f) where xj = j Ax and 
j = 0, ..., J with Ax = 1/ J the uniform grid spacing ( J + 1 grid points). Let the parameter 


Pi = 


a 

a Ax 


(3.38) 


Define the flow residence time in a cell r c = Ax /a (characteristic time due to convection) 
and the time required by the reaction to reach equilibrium r r = 1/a. Then a simple physical 
interpretation of the parameter pi comes from the fact thatpi is equal to the ratio r P /r c . Note 
also that this ratio is the inverse of a Damkohler number. Therefore, the parameter pi is a 
measure of the stiffness of the problem. The smaller pi is, the stiffer the problem becomes. 
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A semi-discrete approximation (for a chosen spatial discretization for the convection and 
source term) of the reaction-convection PDE (3.33a) is then 


IdU 

a dt 


(3.39) 


where U is the vector whose components are u,(<), 1 < j < J. The function F is a 
discrete J -dimensional vector which depends on the grid function U, the parameter p x , and 
the particular spatial finite difference approximation. For simplicity the commonly used spatial 
discretizations (the first-order upwind (UP1), second-order upwind (UP2) and second-order 
central (C2) schemes) are considered for the convection term, and the pointwise evaluation 
(PW), upwind interpolation (UI), and mean average between two neighboring grid points (MA) 
are considered for the (spatial) numerical treatments of the source term. The combination of 
the three numerical treatments of the source term and the three basic schemes used for the 
discretization of the convection term yields 9 spatial finite difference approximations for the 
reaction-convection PDE (3.33a). 


The expressions of the elements fj of F corresponding to the 9 spatial difference approxima- 
tions for (3.33a) and (3.33c) are given below, where we use the obvious notations u-i = u j_i , 
u 0 = and ui = uj+i for the periodic boundary condition (3.33c). 

1. First-order upwind for convection - pointwise evaluation for source term (UP 1PW) 


fj(V) = -p,K - + S{uj). (3.40a) 

2. Second-order upwind for convection - pointwise evaluation for source term (UP2PW) 

fj{U) = -Pi(^ u i - u i-i + + $(«;)• (3.40b) 

3. Second-order central for convection - pointwise evaluation for source term (C2PW) 

f:i u ) = ~^Pi( u j+i ~ «i-i) + S(uj). (3.40c) 

4. First-order upwind for convection - upwind evaluation for source term (UP1UI) 

fj(U) = ~Pi( u j -Uj-i) + 0S(uj-i) + (l-9)S(uj). (3.41a) 

5. Second-order upwind for convection - upwind evaluation for source term (UP2UI) 

fj(U) = “ u j-i + ^ u i-a) + 0S(uj- 1) + (1 - d)S(uj). (3.41b) 

6. Second-order central for convection - upwind evaluation for source term (C2UI) 

fi(U) = -ipi(u i+1 - + 0$(uj_i) + (1 - 9)S(uj). (3.41c) 
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7. First-order upwind for convection - mean average evaluation for source term (UP IMA) 

fi(U) = ~Pi(uj - Uj-i) + ^ S(u j+1 ) + S(uj-i) . (3.42a) 

8. Second-order upwind for convection - mean average evaluation for source term (UP2MA) 

fj(U) = - u i-i + + \ S{uj- 1 ) + 5(u,+i) . (3.42b) 

9. Second-order central for convection - mean average evaluation for source term (C2MA) 

fj( u ) = \ + S{u j+1 ) . (3.42c) 

The parameter 0 associated with the upwind interpolation of the source term in formula (3.41) 
is an extra parameter in the discretization, lying between 0 and 1. For ease of referencing, 
the above 9 spatial discretizations will be denoted by the symbols UP1PW, UP1UI, UP IMA, 
UP2PW, UP2UI, UP2MA, C2PW, C2UI and C2MA. 

Time Discretizations : To contrast the nonlinear behavior between LMMs and explicit Runge- 
Kutta methods, for simplicity, the explicit Euler and the modified Euler schemes are considered. 

Let u” represent an approximation of ti(j Ax, nAf) with a fixed time step Af. Also let U n 
denote a vector with elements u". Then the fully discrete approximation of the PDE (3.33a) 
with explicit Euler time discretization is 

U n+1 =U n + P2F(U n ), (3.43) 

where the parameter pj is simply related to the time step through 

Pi = a At. (3-44) 

With modified Euler time discretization, the fully discrete approximation is 

U n+1 =U n + P2F(U) ; U = U n + ~F(U n ). (3.45) 

In the following the above 18 fully discrete approximations will be denoted for ease of 
reference by the symbols UP1PW/EE, UP1PW/ME, etc. (where EE stands for explicit Euler 
and ME for modified Euler). In all of the computations 0 = c = aj^ in (3.41) is used. 

FDMs Bated on the Characteristic Form : A more physical approach to updating the grid 
value uj at time level n + 1 is to trace back the characteristic passing through (xj, (n 4- 1)A<). 
Denote by (x,nAt) the coordinates of the point on the characteristic at time nAf. Along this 
characteristic, the problem reduces to solving the ODE 
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v r = aS(v) nAt < r < (n 4- 1)A< 


(3.46a) 


with the initial condition 


v(nAt) = u n (i), (3.46b) 

where «”(*) denotes the value of u at point £ and time level nAf, obtained by some 
interpolation method on the adjacent grid function u". For this approach, explicit and implicit 
time discretizations of (3.46a) or, equivalently, explicit and implicit treatments of the source 
term of the PDE (3.33a) are considered. To facilitate the analysis, the convection part of the 
PDE (3.33a) is handled in an explicit way. This means that for the homogeneous part of the 
PDE (3.33a) all the underlying schemes are under the CFL restriction 

a At 

c = ~r— = P 1 P 2 < 1. (3-47) 

An immediate consequence is that £ lies between = (j - 1)A* and Xj = j Ax. Then, 
from a linear interpolation we get 

u” = u n (*) = cu?_, +(1 -cju?. (3.48) 

For fully explicit schemes, the same two explicit time discretizations for the method of lines 
approach are considered here. With explicit Euler, the fully discrete approximation (denoted by 
CHA/EE) takes the form 


u i +1 = «*;-i + (1 - + PiS(u n ). (3.49) 

With the modified Euler time discretization the fully discrete approximation (denoted by 
CHA/ME) takes the form 


«? +1 = + (1 - c)u” + PaS(u), (3.50) 

where 

u = cu"_, + (1 - c)u” + ^5(u n ). (3.51) 

For a less restricted time step, the implicit Euler (IE) and the trapezoidal implicit scheme 
(T) are considered for the source term. With implicit Euler, the fully discrete approximation, 
denoted by CHA/IE, takes the form 

< 3 - 52 ) 

while with the trapezoidal rule, the fully discrete approximation, denoted by CHA/T, takes the 
form 
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“" + ‘ - T 5 (“" +1 ) = cu U + (1 - <0“" + T s ( 4 ”)- 


(3.53) 


The corresponding linearized implicit scheme CHA/LIE associated with the fully discrete 
approximation CHA/IE is given by 


u n+1 = u 1 ? + 
3 3 ' 


CU 7_i -CU? +P35(«7) 


l-p 2 5’(uj) 

with 5' = dS/du = — 2 + 6u — 3u 2 . The linearized trapezoidal method CHA/LT is 

cu - cu? + ^ [5(u?) + 5(« n )] 


(3.54) 


u n+1 = u" + 
3 3 


1 - a S'(«?) 


(3.55) 


3.7.1. Spurious Asymptotes of Full Discretizations 

Besides the three exact steady-state solutions, depending on the numerical methods, either 
the spatial discretizations and/or the time differencing can independently introduce spurious 
asymptotic numerical solutions (see Lafon & Yee (1991) for a detailed analysis). Bifurcation 
diagrams and stability analysis for the exact and spurious asymptotes of the above schemes and 
source term treatments were discussed in Lafon & Yee (1991, 1992). Interested readers are 
referred to these references for more details. 

Since the explicit Euler and the implicit methods are LMMs, no spurious steady state due 
to time discretizations are possible. But, consider for example the various FDMs involving 
modified Euler time differencing (method of lines or characteristic form) and look for the simple 
case of spatially invariant steady states (i.e., the value of uj is the same for all j, 1 < j < J). 
Then it is found that for such FDMs, the value u* of a spatially invariant steady-state must 
satisfy 


u* + ^ S(u *) — u , S(u) — 0. (3.56) 

It can be shown that (3.56) admits the following 9 solutions 
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(3.57) 


in which u* = 0, 1 or 2 are the exact steady-state solutions while the rest are spurious 
steady-state numerical solutions introduced by the modified Euler time discretization. 
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3.7.2. Linearized Behavior vs. Nonlinear Behavior 


To illustrate the differences between the linearized analysis and the nonlinear solution 
behavior. Fig. 3.20 shows the spectral radius around the exact steady- state solution u* = 2 
and the bifurcation diagram obtained with initial data U° = (2., 2.1, 2., 2.1, 2., 2.1, 2.) for the 
scheme UP1UI/EE and p\ = 7. Similar results for the scheme UP1UI/ME are shown in Fig. 
3.21. For pi = 7, the scheme UP1UI/ME exhibited two disjoint linearized stability intervals. 
From Fig. 3.21 we observe that outside these stability intervals the scheme does not necessarily 
diverge (as indicated by the linearized analysis) but can converge to a spurious asymptotic 
solution. 

Another example which shows the importance of nonlinear analysis is that of two schemes 
that exhibit the same linearized behavior yet have different nonlinear behavior (true behavior). 
For example, the linearized stability analyses for schemes UP1UI/EE and CHA/EE (see Figs. 
5a and 7a of Lafon & Yee (1992)) are identical even though the bifurcation diagrams obtained 
with the same initial data and pi = 7 are different (see Figs 19a and 20a of Lafon & Yee 
(1992)). 


3.7.3. Spurious Steady States and Nonphysical Wave Speeds 

The possible connection of the numerical phenomenon of incorrect propagation speeds of 
discontinuities with the existence of some stable spurious steady states introduced by the spatial 
discretization is discussed here. The boundary condition (3.33c) and the following piecewise 
constant initial data 


u 


o 



2 

0 


0 < X < Zrf, 

Xj < X < 1 


(3.58) 


are considered. The constant value xj denotes the location of the discontinuity. The exact 
solution of (3.33a,b,c) with initial data (3.58) is simply 


u(*,<) = u°(* — af) (3.59) 

which is a wave traveling at constant speed a. 

For an explanation of how numerical methods applied to this model PDE are likely to produce 
wrong speeds of propagation for the initial data (3.58), the reader is referred to LeVeque & 
Yee (1990). This phenomenon lies with the smearing of the discontinuity caused by the spatial 
discretization of the advection term. This introduces a nonequilibrium state into the calculation. 
Thus, as soon as a nonequilibrium value is introduced in this manner, the source term turns on 
and immediately tries to restore equilibrium, while at the same time shifting the discontinuity 
to a cell boundary. 

For simplicity, consider the first-order upwind spatial discretization with the explicit Euler 
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time discretization for (3.33a, c) and (3.58). Assume an equal spatial increment of J intervals so 
that the discretized initial data associated with (3.59) is 


| 2 1 <j<K 


(3.60) 


with the index K depending on the constant z*. In addition, assume that the convection speed 
a = 1 so that At = c/J. In the computation J = 20 and the Courant number c - .05. With 
these assumptions, only the stiffness a of the source term is a free parameter. Define the average 
wave speed W for the numerical solution as follows 


W = 


Ax 

2nAt 


r J+i 

L j=i J 




(3.61) 


where the average is taken on the time interval 0 to nAt. Figure 3.22 shows this average speed 
versus pi (proportional to 1/a). It reveals the fact that when pi becomes large (or, similarly 
when the source term is not stiff) the numerical wave speed tends to the exact solution (a = 1), 
while for pi < .25 the numerical solution is a standing wave (the average speed being 0). 

This zero wave speed is indeed a stable steady-state solution of the discretized equation 
but not a solution of the continuum PDE. Since the explicit Euler time differencing has 
the property of not producing spurious steady-state numerical solutions, this spurious stable 
steady-state numerical solution must have been introduced by the spatial discretization. In 
fact, it is evident from the bifurcation diagram shown in Fig. 2 of Lafon & Yee (1991) 
that this spurious steady state lies on the stable branch originating at pi = 0 from the state 

tij = 2, ...., uj, = 2, tijc+i = 0, ...., u j - 0. 


3.7.4. Numerical Basins of Attraction 


In order to show the global nonlinear solution behavior of the schemes, numerical basins 
of attraction (of the underlying schemes) are compared with each other as well as with the 
exact basin of attraction for u* = 2 (BA( 2)). Due to the complexity and CPU intensive nature 
of the computation, unlike the detailed basins of attraction presented in Yee & Sweby (1994, 
1995a,b), only a fraction of the basins of attraction are computed. 

To compute a partial view of some numerical basins of attraction for u* = 2 (BA( 2)), a set 
of initial data in a two-dimensional plane was selected. Even with this restriction, the analysis 
is still very complex and the computation is CPU intensive; it was performed only for the case 
J = 4 (5 grid points). The equation of the chosen plane is 

Ul = 2 ; u 4 = 2. (3.62) 

A set of 121 initial data in the plane (3.62) were obtained by confining u 2 and in the square 
1.1 < u 2 < 3.1, 1.1 < uj < 3.1 with a uniform increment of At* 2 = Attj = 0.2. For each 
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initial datum, 5000 preiterations were performed. The asymptotic behavior is plotted in the 
(u 2 , uj) plane. In all of the figures, open triangles belong to the numerical basin of attraction 
BA(2). Dark squares belong to the numerical basin of attraction of a spurious steady state. 
Dark circles are the numerical basins of attraction of other (exact or spurious) asymptotic 
solutions, while a blank space denotes a divergent solution. Note that the whole square domain 
1.1 < u 2 < 3.1, 1.1 < uj < 3.1 is inside the true basin of attraction of the exact steady 
state u* = 2 and therefore the true behavior should be an open triangle for all the initial data 
considered. The study in Lafon & Yee ( 1991) indicated that the modified Euler time differencing 
has a smaller attractive region than the explicit Euler for the domain considered even though 
in terms of linearized stability and accuracy, the modified Euler exhibits an advantage over the 
explicit Euler. 

For example, for pi = 0.1 and pj = 0.5, implicit treatments of the source term exhibit a 
larger attractive region of initial data than the explicit treatments of the source term. However, 
as the time step is further increased, an adverse behavior is observed contrary to the common 
belief that implicit schemes can operate with much higher time step and still produce the 
correct steady-state numerical solution. Since the source term is handled implicitly, the classical 
guideline for the time step constraint is given by the explicit discretization of the convective term 
(homogeneous part of the PDE (3.33a)), or equivalently by the CFL constraint c = pipi < 1. 
Thus, the permissible time step for the implicit treatments of the source term is larger than 
explicit treatments of the source term. However, it is evident from the computation shown on 
Figs. 3.23 - 3.25 for implicit schemes CHA/LIE and CHA/LT with a Courant number set equal 
to 0.3 (pa = 3) and 1 (pj = 10) respectively, that these implicit schemes no longer give the 
correct asymptotic solution, in particular for the scheme CHA/LT. 

3.8. Spurious Dynamics in Time-Accurate (Transient) Computations 

In the examples chosen by Lorenz (1989), he showed that numerical chaos always precedes 
divergence of a computational scheme. He suggested that computational chaos is a prelude 
to computational instability. Poliashenko & Aidun (1995) showed that this is not a universal 
scenario. In previous sections, we have shown that numerics can introduce chaos. Using a simple 
example, Corless (1994b) showed that numerics can suppress chaotic solutions. The work of 
Poliashenko and Aidun also discussed spurious numerics in transient computations. Adams et 
al. (1993) discussed spurious chaotic phenomena in astrophysics and celestial mechanics. They 
showed that the source of certain observed chaotic numerical solutions might be attributed to 
round-off errors. Adams also discussed the use of interval arithmetics (interval mathematics or 
enclosure methods) to avoid this type of spurious behavior. Moore et al. (1990) discussed the 
reliability of numerical experiments in thermosolutal convection. Keener (1987) discussed the 
uses and abuses of numerical methods in cardiology. 

It is a common misconception that inaccuracy in long-time behavior of numerical schemes 
poses no consequences for transient or time-accurate solutions. This is not the case when one 
is dealing with “genuinely nonlinear DEs (Jackson 1 989)’ ’. For genuinely nonlinear problems. 
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due to the possible existence of spurious solutions, larger numerical errors can be introduced by 
the numerical methods than one can expect from a local linearized analysis or weakly nonlinear 
behavior. Lafon& Yee (1991) illustrated this connection. Section 3.7.3 summarizes their result. 
The situation can become more intensified if the initial data of the DE is in the basin of attraction 
of a chaotic transient (Grebogi et al. 1983) of the discretized counterpart. In fact, it is possible 
that part or all of the solution trajectory is erroneous. Section V shows a practical example of a 
chaotic transient near the onset of turbulence in direct numerical simulations of channel flow by 
Keefe (1996). Since numerics can introduce and suppress chaos, and can also introduce chaotic 
transients, the dangers of relying on numerical tests for the onset of turbulence and chaos is 
evident. 


There has been much debate on the overall accuracy away from shocks of high- resolution 
shock-capturing methods. Donat (1994) addressed this issue from a theoretical standpoint while 
Casper & Carpenter ( 1 995) illustrated it with a shock induced sound wave model. Casper and 
Carpenter concluded that there is only first-order accuracy downstream of the sound-shock 
interaction using a spatially 4th-order ENO scheme. Sections 5.2 and 5.3 illustrate two additional 
types of spurious numerics in transient computations. 


IV. Spurious Dynamics in Steady-State Computations 


Any CFD practioner would agree that making a time-marching CFD computer code 
converge efficiently to a correct steady state for a new physical problem, that was previously 
not understood, is still an art rather than a science. One usually has to tune the code and 
still encounters one or more problems such as blow-up, nonconvergence, nonphysical, or 
slow convergence of the numerical solution. Some of these phenomena have been reported 
in conference proceedings and reference journals, but the majority have been left unreported. 
Although these behaviors might be caused by factors such as poor grid quality, under-resolved 
grids, improper numerical boundary conditions, etc., most often they can be overcome by 
employing standard procedures such as using physical guidelines, grid refinement, improved 
numerical boundary treatments, halving of the time step, and using more than one scheme 
to double check if the numerical solution is accurate and physically correct. However, 
these standard practices alone may sometimes be misleading, not possible (e.g., too CPU 
intensive) or inconclusive due to the various numerical uncertainties (see section I) that can be 
attributed to the overall solution process. Consequently, isolation of the sources of numerical 
uncertainties is of fundamental importance. Section III isolates some of the spurious numerics 
for elementary models. Complementing the phenomena observed in Section III, this section 
illustrates examples from CFD computations. These examples were published in Yee & Sweby 
(1996). We concentrate mainly on the convergence issues that are contributed by the spurious 
dynamics that are inherent in the schemes. Section 4.2. 1 was written by Bjorn Sjogreen of the 
Uppsala University, Sweden. 


47 



4.1. A 1-D Chemically Relaxed Nonequilibrium Flow Model 4 

The model considered is a nonequilibrium flow field relaxation of the one-dimensional Euler 
equations fora (N 2 , N) mixture. This type of flow is encountered in various physical situations, 
such as shock tube experiments (the mixture behind the shock being in a highly nonequilibrium 
state) or an hypersonic wind tunnel. Under these assumptions the model can be expressed as a 
single ODE, 

^ = S(p,T,z), (4.1) 

where z is the mass fraction of the N 2 species, p is the density of the mixture and T is the 
temperature. There are two algebraic equations for p and T. 

This system consists of a disparity range of parameter values (many orders of magnitude) 
and is stiff and highly nonlinear. Even when care is taken, the following features are observed 
when Runge-Kutta schemes are used for the integration (Sweby et al. 1995): 

• restricted basins of attraction 

• spurious equilibrium values 

• oscillatory solutions 

The derivation of the model is as follows. The one-dimensional steady Euler equations for a 
reacting (N 2 ,N) mixture are 


s («*>“) = 

(4.2.) 

frl*. 

II 

O 

(4.2b) 

O 

II 

+ 

■*1-8 

(4.2c) 

d r i 

dx + p) = o, 

(4.2d) 


where (4.2a) is the balance equation for the N 2 species and wpi a is the production rate of the 
N 2 species with density pn 2 . Equation (4.2b) is the continuity equation, (4.2c) the momentum 
equation and (4.2d) the energy equation, in which p, u, E and p are density, velocity, total 
internal energy per unit volume, and pressure, respectively. 

The production rate w jv, of species N 2 is the sum of the production rates for the two reactions 


4 We would like to acknowledge Andre Lafon for the original formulation and the earlier study used 
in this section; presented at the ICFD Conference on Numerical Methods for Fluid Dynamics, April 
3-6, 1065, Oxford, UK. 
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(4.3a) 

(4.3b) 


N 2 + N 2 2N + N 2 
N 2 + N 2N + N 

and are computed using the Park’s model (Park 1985) that has been largely used for hypersonic 
computations. See the Workshop on Hypersonics (1991) for some discussion. These reaction 
rates involve an equilibrium constant, K eq (see below), which is determined by a polynomial 
fitting to experimental data, and as such is only valid for a certain range of temperatures. In 
particular, a cut-off value has to be introduced for low temperatures, a typical choice being 
Tmin = lOOOif (Mulard & Moules 1991). 

The system (4.2) and (4.3) must be closed by a thermodynamical representation of the 
mixture. Here a simple model, with no vibrational effects, has been chosen. The details have 
been omitted for brevity. 

Equations (4.2b)-(4.2d) simply integrate to give 


pu = qoc, (4.4a) 

pu 2 + P = P<x>, (4.4b) 

H = *L±P = J 0#) (4.4c) 

P 

where H is the total enthalpy and qoo, Poo and Hoo are all constants. Finally, denoting the mass 
fraction of the N 2 species by 


z = 


PNj 

P 


(4.5) 


and using the Park’s reaction rate model and the thermodynamical closure, (4.1) can be written 
as 


dz 

dz 


= S(p,T,z) 

x [ aAiz(l - z) 2 — .Aiz 2 + 2aA 2 (l — z) s — 2A 2 z(l — z)] , 


(4.6a) 


where 


a = 


M\K t 




10 4 

K tq = 10* exp [cj + c 2 Z + cjZ 2 + c^Z 1 + C5Z*], Z = ~jT- (4.6b) 
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The density p is obtained from 


*L(8 - 2 *)Q) -(10-3z)Poo^ +2(2-z)[H 00 -(l-z)el] =0 (4.6c) 

and the temperature T from 


with pressure 


The model uses the constants 


' Cl = 3.898 ■ 


Ai = 3.7 x 10 15 ‘ 

c 2 = -12.611 


A 3 = 1.11 x 10 ia 

cj = 0.683 


B - -1.6 

c 4 = -0.118 


0 = 1.132 x 10 5 

c 5 = 0.006 


e\ = 3.355 x 10 7 

.Mi = 28 x 10* . 


. R = 8.3143 . 


where qoo, P and Hoo are input parameters set to be 0.0561, 158,000 and 27,400,000, 
respectively, in our calculations. A limitation of the model is T > TL;„ = lOOOJif. In addition, 
solutions are nonphysical if z £ [0, 1], if p < 0 or if p is complex. 

Equation (4.6a) was integrated using the Euler, modified Euler, improved Euler, Heun (R-K 
3), Kutta (R-K 3) and R-K 4 Runge-Kutta schemes. In each case the computations were 
performed for a range of initial z and integration steps A*. The discretized equations were 
preiterated 1000 steps before a full bifurcation diagram together with basins of attraction were 
produced. (The computations were carried out on the CM- 200 Connection Machine at the 
Edinburgh Parallel Computing Centre, UK.) 

There are two strategies possible when implementing the Runge-Kutta schemes. One is to 
freeze the values of p and T at the beginning of each step when calculating S(p, T, z) at the 
intermediate stages. The other is to update the values at each evaluation of the function S. The 
results presented here employ the latter strategy since this is the more proper implementation; 
however, it is interesting to note (see below) that results obtained by freezing p and T for 
intermediate calculations exhibit a slightly richer dynamical structure. 

Figure (4.1) shows the results obtained from our calculations. In all of these plots the shaded 
region denotes the basin of attraction in which combinations of initial z values and step size 
Ax converge to the fixed point, depicted by the solid black line. The unshaded regions indicate 


T = 


Mip 


R ( 2 - z)p 


__ p _ 
P *oo 


(4.6d) 


(4.6e) 
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regions where the combinations of z and A* do not converge to a physical solution of the 
problem. 

As can be seen in all cases the basins of attraction narrow considerably for the larger values 
of A*. (Note that the axis scale is 10 — 5 !) The explicit Euler scheme (Fig. 4.1a) obtains 
the correct fixed point up to its stability limit, where there is a very small region of period 
2 (oscillatory) behavior before it becomes nonphysical. Similar behavior is observed for the 
improved Euler (Fig. 4.1c) and Kutta (Fig. 4.1e) schemes, the latter also exhibiting a much 
more constricted basin of attraction for any given A*. The Heun scheme (Fig. 4. Id) exhibits a 
distinct region of oscillatory behavior just above the linearized stability limit. 

As is typical with the modified Euler scheme (Fig. 4.1b), a transcritical bifurcation occurs 
at its stability limit which leads to a spurious (A* dependent) solution. Note also the solid 
line at about z = 0.25 down on the plot, outside of the shaded region. This appears to be an 
unstable feature picked up by our method of fixed point detection (comparison of initial data 
with the 1000th iterate) and is unlikely to arise in practical calculations unless the initial data 
are on this curve. The R-K 4 scheme (Fig. 4. If) also exhibits a transcritical bifurcation at the 
linearized stability limit; however this is discernible more by the sudden narrowing of the basin 
of attraction since the spurious fixed point varies only slightly with Ax. 

Finally we note that if the values of p and T are frozen for intermediate calculations the 
dynamics are somewhat modified. All schemes with the exception of explicit Euler have a 
slightly larger basin of attraction for values of Ax within the stability limit and all schemes 
have period two behavior at the stability limit, there being no transcritical bifurcations for any 
of the schemes. The modified Euler scheme also has embedded period doubling and chaotic 
behavior below the linearized stability limit. 

What we have shown with this example is that, although the dynamical behavior is perhaps 
not as rich as in some of our simple examples, spurious features can still occur in practical 
calculations and so care must still be taken in both computation and in interpretation. 

4.2. Convergence Rate & Spurious Dynamics of High-Resolution Shock-Capturing Schemes 

We have seen in Section III elementary examples and references cited therein on how the 
proper choice of initial data and the step size combination can avoid spurious dynamics. Yet 
for other combinations the numerical solutions can get trapped in a spurious limit cycle. We 
have also seen that the convergence rates of the schemes are greatly affected by the step sizes 
that are near bifurcation points. Here we include the dynamics of full discretization of two 
nonlinear PDE examples. The spatial discretizations are of the high-resolution shock-capturing 
type (nonlinear schemes). This includes TVD and ENO schemes. Section 4.2.1 discusses how 
this nonlinear scheme affects the convergence rate of systems of hyperbolic conservation laws. 
Section 4.2.2 illustrates the existence of spurious asymptotes due to the various flux limiters 
that are built in to TVD schemes. 
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4.2.1. Convergence Rate for Systems of Hyperbolic Conservation Laws 


This section summarizes the results of Engquist & Sjogreen (1995) and Sjogreen (1996, 
private communication). These results concern the convergence rate for discontinuous solutions 
of a system of hyperbolic conservation laws. For a scalar nonlinear conservation law, the 
characteristics point into the shock. According to the linear theory of Kreiss & Lundqvist 
(1968), dissipative schemes damp out errors propagating backwards against the direction of the 
characteristics. Thus it is reasonable to expect that the locally large errors at the shock stay in 
a layer near the shock. In numerical experiments we usually obtain 0(h p ) convergence away 
from the shock with difference schemes of formally p th order. 

For the systems case, the same reasoning from the scalar conservation law cannot be applied. 
In this case, we have other families of characteristics intersecting the shock causing the situation 
to be more involved. Thus it is possible that the large error near the shock propagates out into 
the entire post shock region by following a characteristic which emerges from the shock. 

This effect cannot be seen in a simple scalar Riemann problem (problem with jump initial 
data), because exact global conservation determines the post shock states. The system model 
problem, taken from Engquist & Sjogreen (1995), 


u t + (u a /4) m = 0, -oo < * < oo, 0<< (4.7a) 

v t +v m +flr(u) = 0, s(u) = (u + l)(u - l)(l/2 - u) (4.7b) 

gives an example of propagation of large errors. The function g(u) has the properties 
0(1) = ff(— 1) = s(l/2) = 0, and g(u) ^ 0 for -1 < u < 1 except at u = 1/2. The initial 
data was given as 


, x r 1 * < 0 

“•(•) = | _! x>Q > «•(*) = 1 (4.7c) 

so that the exact solution of the u equation is a steady shock. The eigenvalues of the Jacobian 
matrix of the flux (u 2 /4, v) T for (4.7) are Ax = u/2 and X 2 = 1. The eigenvalue Ax = u/2 
corresponds to a strictly nonlinear field, and A 2 = 1 corresponds to a linearly degenerate field. 

With this initial data (4.7c), it gives rise to a steady 1-shock, with the 1 -characteristics having 
a slope 1/2 to the left of the shock and a slope - 1/2 to the right of the shock. These thus intersect 
the shock when time increases. The 2-characteristics of the linear field, have slope 1 on both 
sides of the shock. These characteristics thus enter the shock from the left and exit to the right. 
The v component of the solution is passively advected along the 2-characteristics. When these 
characteristics exit from the shock at * = 0, an error, coming from poor accuracy locally at the 
shock, is picked up and advected along with the solution into the domain * > 0. The shock 
curve * = 0 ( in the x-t plane ) acts as an inflow boundary for the domain x > 0. The error 
coming from the shock is similar to an error in given inflow data, and is therefore not affected 
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by the numerical method used in the interior of the domain. Thus is not surprising that this 
error is of first order, even when the equation is solved by a method of higher formal order of 
accuracy. 

Figure 4.2 shows the numerical solution, computed by a second-order accurate ENO method 
using 50 grid points at the time T = 5.68. The points in the shock give a large error which is 
coupled to the v equation through g(u). The exact solution for v is 1. Numerical investigation 
of the convergence rate of the error in v to the right gave the exponent 1 .047. One thus has 
first-order convergence for this second-order accurate method. 

Similar effects can be seen in computing the quasi one-dimensional nozzle flow. Engquist 
& Sjogreen (1996) computed the solution on the domain 0 < * < 10 for a nozzle with the 
following cross sectional area variation 

A( x) = 1.398 + 0.347 tanh(0.8* - 4). (4.8) 

This problem is studied in Yee et al. (1985) for a class of explicit and implicit TVD schemes. 
The solution has a steady shock in the middle of the domain. Figure 4.3 shows the error in 
momentum for the steady-state solution on grids of 50, 100, and 200 points for a fourth-order 
ENO scheme and a second-order TVD scheme. For the fourth-order method, the convergence 
exponent is 3.9 before the shock and 1.0 after the shock, when going from 100 to 200 points. 
For the second-order TVD the same quantities have the values 2.2 and 1.1 respectively. 

Sjogreen (1996) recently conducted the same numerical study for the two-dimensional 
compressible Euler equations for a supersonic flow past a disk with Mach number 3. The 
equations were discretized by a second-order accurate uniformly nonoscillatory (UNO) scheme 
(Harten 1986), which unlike TVD schemes, is formally second-order everywhere including 
smooth extrema. He computes the error in entropy along the stagnation line for the steady-state 
solution on grids with 33 x 17, 65 x33, and 129x65 grid points. The result is shown in Fig. 4.4, 
where the error and convergence exponent in the region behind the bow shock are plotted. The 
convergence exponent is between the 65 x 33 and 129 x 65 grids. The disk has radius 0.5, and 
it is centered at the origin, which means that the line is attached to the wall for —0.5 < x. A 
convergence exponent of 1 .5 is observed for this formally second-order method. 

4.2.2. Spurious Dynamics of TVD Schemes for the Embid et al. Problem 5 

It has long been observed that the occurrence of residual plateauing is common when TVD 
and ENO types of schemes are used to time march to the steady state. That is, the initial 
decrease in the residual levels out and never reaches the convergence tolerance. See Yee (1986, 
1989) and Yee et al. (1990) for some discussion. This has often been overcome by ad hoc 
modification of the flux limiter or similar device in problem regions. 


1 We would like to acknowledge Paul Burton for the computations used in this section. 
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A recent study (Burton & Sweby 1995) investigated this phenomenon using a dynamical 
systems approach for the one-dimensional scalar test problem of Embid et al. (1984) 

«t+Q« 2 ) = < 7 (*)u = (6x - 3)u, *<=(0,1), (4.9a) 

with boundary conditions 

u(0) = 1, u(l) = -0.1. (4.9b) 

This equation with the flux function f(u) = u a /2 has the property that there are two entropy 
satisfying steady solutions, consisting of stationary shocks jumping between the two solution 
branches 


u l (x) = Zx(x - 1) + 1 (4.10a) 

u r (*) = 3x(* — 1) — 0.1. (4.10b) 

For this problem the two possible solutions consist of a single shock, either approximately at 
*i = 0.18 or xj = 0.82. It can then be shown (see Embid et al. 1984 for details) that the 
solution with a shock at x\ is stable to perturbations while the solution with a shock at x 2 is 
unstable. 

Embid et al. solved (4.9) using three different methods, the first-order implicit upwind 
scheme of Engquist and Osher, its second-order counterpart, and the second-order explicit 
MacCormack scheme. All three schemes used time stepping as a relaxation technique for 
solving the steady- state equation. The initial conditions were taken to follow the solution 
branches (4.10) from the boundary values, with a single jump between the two branches. The 
results obtained showed that, although the implicit schemes allowed large time steps and hence 
fast convergence, if the initial jump was taken too near the unstable shock position x 2 , then for 
some ranges of Courant number, 


c = 



(4.11) 


the schemes would converge to the physically unstable shock. This phenomenon was studied 
both for these three schemes and a variety of flux limited TVD schemes (Sweby 1984) in Burton 
& Sweby (1995), where not only the full problem was studied but also a reduced 2x2 system 
was investigated using a dynamical system approach. We report on this investigation here. 


The schemes investigated were explicit and implicit versions of the Engquist-Osher and 
TVD flux limiter schemes using the minmod, van Leer, van Albada and superbee flux limiters. 
In all cases a method of lines implementation was adopted, the spatial discretization being 
performed separately from the time discretization. (This is necessary so that the final state 
does not depend on the Courant number, which would be the case if a Lax-Wendroff type 
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discretization were used.) For the time discretization, forward Euler was used for the explicit 
implementations while linearized implicit Euler was used for the implicit computations. For 
the second-order flux limiter schemes the Jacobian matrix used was taken to be that of the 
first-order Engquist-Osher in order to allow easy inversion. 

The schemes are 

(a) Explicit Scheme 


“" +1 = «” - AA_(/r +1 + /+) + A (,(«)«? 

- iAA_(<Hr+)(A/ i+ })+ -*(rT +I )(A/ j+| )-], (4.12) 

(b) Linearized Implicit Scheme 


J(u?)(«; +1 - «J] = - AA_(/r +1 + ff) + 


- -AA_(«r;)(A/ J+ .)+ -^r7 +1 )(A/ i+ ,)-], 


where ff are the Engquist-Osher numerical fluxes 


ff = /( mar(ttj, 0)), 
fj = /(min(u J -, 0)). 


The flux differences are given by 


(A/ i+ i ) + = - (ff +t - /(u i+1 )), (A/ i+ i) = (f j - /(«;)), 

and the solution monitors by 




Wi-i ) 


±1 


UA/z+i)* J • 

Finally, J is the Jacobian matrix and the flux limiter <j>(r) is one of 


(4.13) 

(4.14a) 

(4.14b) 

(4.15) 

(4.16) 


<f>o(r) — 0 — first-order Engquist-Osher (E-O) scheme, 

= max(0,min(r, 1)) — the minmod limiter, 

— max(0, min(2r, 1), min(r, 2)) — Roe’s superbee limiter. 


<t>VL(r) 
<t>VA{r ) 


van Leer’s limiter, 

l + |r 


— van Albada’s limiter. 

1 + r 2 


(4.17) 

(4.18) 

(4.19) 

(4.20) 

(4.21) 
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The experiments reported in Burton & Sweby (1995) used a grid spacing of Ax = 0.025 
with initial conditions consisting of a single jump between the solution branches (4.10) near 
either the stable shock (x t ) or the unstable shock (x 2 ). The convergence criterion used was the 
following bound on the residual 


-u?| < IQ' 15 , (4.22) 


with an upper limit of 2000 iterations being performed. 

The results of applying these schemes to problem (4.9) largely echoed those reported by 
Embid et al. For the explicit schemes convergence, when it occurred, was to the stable shock. 
It was found that there were regions of Courant number (c < c<) for which the schemes 
converged, regions (c» < c < cj) for which convergence did not take place using (4.22) within 
2000 iterations, and regions (c > cj) for which the schemes were unstable. This is summarized 
in Table 4. 1 . The absence of an entry in Table 4. 1 corresponds to residual plateauing. 

Notice that for the superbee flux limiter there was no range of Courant numbers for which the 
scheme converged. Closer inspection reveals residual (defined as r n = |u n+1 — u n |) plateauing 
at around 10 -3 . For the other schemes when c,- < c < Cj the nonconvergence observed arises 
from a similar process, except that the residual does not necessarily level out completely, but 
decreases at a very gradual rate, resulting in very slow convergence. 

The implicit scheme experiments revealed that the choice of initial conditions could cause 
convergence to the unstable shock for certain ranges of Courant number. For an initial jump 
near the stable shock, the schemes (with the exception of the van Leer limiter) converge to 
the stable shock for c < 11. However, for an initial discontinuity near the unstable shock, 
convergence could sometimes be towards the unstable shock. The situation is summarized in 
Table 4.2, where again the absence of an entry corresponds to residual plateauing. 

To gain further insight into this problem. Burton & Sweby (1995) considered a reduced 
problem consisting of two free points at one of the shocks, with exact solution values being 
imposed as boundary conditions on either side. This then leads to a two dimensional dynamical 
system which, although obviously a gross simplification of the full problem, was hoped to still 
maintain some of the qualitative behavior. 

The situation is as shown in Fig. 4.5, where the free points are X and Y, the remaining 
points (Uu, Ui , U rr and U r ) being set at exact analytic values to provide boundary conditions. 
Two such values are needed on either side to provide the necessary information for the flux 
limiters. Substitution of these points into the numerical scheme then leads to a two-dimensional 
system. For example, the explicit Engquist-Osher scheme yields 
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x n+1 = x n - — 

40 


yn-fl y-n 


At 

40 


f~(Y n ) + /+(*”) - /“(X n ) ~ f + (Ui) 
r(Ur) + f + {Y n ) - f-(Y n ) - f + (X n ) 


39 

— A<X n 
20 


36 

— A<y n , 

20 


(4.23a) 

(4.23b) 


where a step size of A* = ^ has been used. 

For the first-order explicit and implicit schemes some analysis on the reduced problem can be 
performed. Table 4.3 summarizes the findings. Note that for both schemes spurious fixed points 
are introduced by the simplification of the problem. These both have X and Y of the same 
sign and would not be tolerated for the full problem. It is only the proximity of the boundary 
conditions for the reduced problem which allow them to exist as fixed points. However, the 
remaining fixed points and their stability agree well with numerical results obtained for the full 
problem. 

Analytical results could only be obtained for the first-order scheme and so numerical 
experiments were performed for the flux limiter schemes. These consisted of generating 
bifurcation diagrams for X and Y against At and the plotting of basins of attraction in the 
(X, Y) plane for fixed values of At. The explicit schemes were shown to possess no spurious 
dynamics below their respective stability limits, apart from that introduced by the simplification 
of the problem (i.e. outside of the quadrant (X > 0, Y < 0)). As At was increased above the 
stability limit the schemes entered a period of bifurcation and chaos accompanied by a dramatic 
shrinkage in the numerical basins of attraction. 

The dynamics of the implicit schemes at the unstable shock showed the falsely stable fixed 
point becoming stable for large values of At. For all the limiters tested the stabilizing of 
the fixed point was accompanied by the introduction of additional, spurious (period 2) fixed 
points. These spurious solutions caused a reduction in size of the basin of the falsely stable 
fixed point. The fact that the more compressive limiters took longer to recover from the effects 
of the spurious fixed points seems a possible cause of the phenomenon of residual plateauing 
experienced in the full problem. Due to a space limitation, see Burton & Sweby (1995) for the 
illustrations. 

It must be realized that although the residual plateauing illustrated is around the physically 
unstable shock (to which we would usually not wish to converge), the fact that it is not a 
repelling phenomenon will in itself have repercussions on convergence to the correct, physically 
stable shock. Indeed no such nonconvergence was observed for Courant numbers greater than 
1 1 . We conclude this section by emphasizing that the reduced problem indicated a possible 
cause for residual plateauing. However, for certain situations the dynamics of the full problem 
does not coincide precisely with that of the reduced problem. 
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4.2.3. The Dynamics of Grid Adaption 


Consider a model convection-diffusion equation of the form 

«t + /(«)» = eu„. (4.24) 

with the linear case, /(u) = u and the nonlinear case, /(u) = |u 2 (the Burgers equation) 
respectively. The boundary conditions for the linear case are u(0,<) = 0 and u(l,<) = 1 
and for the nonlinear case, u(0,<) = 1 and t*(l,<) = -1. These boundary conditions result 
in steady-state solutions of a boundary layer at * = 1 and a viscous shock at z — 1/2, 
respectively. In both cases the steepness of the feature is governed by the diffusion coefficient 
parameter e. Besides its steepness feature, one of the main reasons for considering the linear 
convection-diffusion equation is to show that grid adaptation alone and/or nonlinear schemes 
such as TVD schemes can introduce unwanted dynamics to the overall solution procedure. The 
authors realize that model (4.24) is not the best model to illustrate the dynamics of the studied 
schemes since the model is not stable under perturbations. However, it serves to show what 
type of spurious numerics would occur under such an environment. 

One common criterion used for grid adaptation is the equidistribution of a positive definite 
weight function t »(*,<), often taken to be some monitor of the numerical solution u(x,<) of 
the underlying PDE. A grid *«< *i (<)<••• < < xj, where x 0 and xj are fixed, 

equidistributes w(x, <) (at time t) if 

f m i 2 f m J 

I w(x,t)dx = I w(x,t)dx = — I w(x,t)dx , (4-25) 

••i - 1 *•/ * *»o 

for j = 1, . . . , J. A one-parameter family of weight functions 

w(x, t) = y/l — a + au*(*,<), a E [0,1], (4.20) 

can be chosen where a = 1/2 corresponds to equidistribution in the arc-length and a = 0 
yields a uniform grid. Approximating w(x, t) to be constant in each interval (zj-i ,*>) yields 


(v/l-a + awj) a (*i - *j- 1 ) = (yi-a + aul) * (x j+1 - *_,). (4.27) 

Given a numerical solution of the PDE we can approximate the derivatives by 

u.li- 1/2 » ~ Ui ~ 1 . (4.28) 

Xj — Xj-i 

Equation (4.27) is nonlinear in {xj} if we use (4.28). However, (4.27) is linear in {xj} if {xj} 
in (4.28) uses the existing grid. In this case we can solve the tridiagonal system (4.27) for a new 
set of {xj} to obtain an updated grid. 
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Given a set of initial data and an initial grid, the procedure is to numerically solve the PDE 
and (4.27) in a time-lagged manner. We use nodal placement and the I 2 norm of the solution 
to illustrate our results. We use the previous time step value for xj in (4.28) to achieve a linear 
tridiagonal system for the updated grid in (4.27). Our preliminary study shows that the solution 
procedure of Ren and Russel (1992), and Budd et al. (1995a) in solving the coupled PDE and 
(4.27) are less stable than the present linearized form. See also Neil (1994) for a similar study 
and conclusion. The regridding strategy adopted was to regrid after every time step of the 
PDE method, either interpolating updated solution values from the old grid or performing no 
adjustment at all due to grid movement. This latter approach in effect presents the PDE method 
with new initial data to the problem at each step. 

The dynamics of the above one-parameter family of mesh equidistribution schemes coupled 
with different spatial and time discretization were studied numerically in Sweby & Yee (1994) 
and Yee & Sweby (1995a,b) using the above numerical procedure. The spatial discretizations 
include three-point centred, second-order upwind and second-order TVD schemes. The time 
discretizations include explicit Euler, second- and fourth-order Runge-Kutta and the linearized 
implicit Euler methods. In a parallel study, Budd et al . ( 1 995b) made use of the AUTO computer 
bifurcation package (Doedel 1986) to obtain bifurcation diagrams for similar grid adaptation 
methods for the steady part of the above PDEs with a different form than (4.27). However, 
the dependence on known solutions of the discretized PDEs and grid equations as starting 
values limits its usage. In Sweby & Yee (1994) we utilize the power of the highly parallel 
Connection Machine CM-5 to undertake a purely numerical investigation into the dynamics of 
the time-marching adaptive procedure. 

We divided a chosen parameter space (e.g., e) into 512 equal intervals, with all other 
parameters (a, At, initial data) fixed. For each chosen parameter value, we iterated the 
discretized PDE and the grid function, in general, 4000 steps (8,000 steps for explicit methods) 
to allow the solution to settle to an asymptotic state. Then we performed a series of time 
step/regridding stages, during which we investigated the dynamics by producing an overlaid 
plot of the I 2 norms of the numerical solution and the grid distribution at each step. This 
resulted in a bifurcation type diagram or the grid displacement diagram as a function of the 
physical or discretized parameters. We also performed numerical studies by only preiterating 
the discretized PDEs to the steady state for a fixed grid before solving both the discretized PDE 
and the grid adaptation function. We found in most cases the solution process is less stable and 
more likely to get trapped in a spurious mode than in the aforementioned process. 

For this study, we took into consideration the grid density, an even and odd number of nodes, 
and whether or not there is interpolation after each regridding. Tlie grid density studies consist 
of 4, 5, 6, 9, 10, 19, 20, 49 and 50 grid points. There is no apparent sign of even or odd grid 
dependence. The resolution and stability of TVD schemes are also grid independent. However, 
the central difference scheme experienced instability more often for coarser grids, and the 
second-order upwind is slightly more stable, with better resolution than the central scheme. As 
expected, the stable time step required for the explicit methods was orders of magnitude lower 
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than that for the implicit method. For the TVD schemes, comparison of the dynamical behavior 
of the five limiters of (3.50) of Yee (1989) was performed. Four of the limiters are the same as 
(4.17)-(4.21). Due to the simplicity of the PDEs, their dynamical behavior is similar, although 
there were slight differences in the stability and resolution. Due to space limitations, we 
summarize the results without presenting the actual computational figures. Interested readers 
should refer to our original papers for details. 

Summary : We consider separately the cases “with” and “without” interpolations. The term 
“scheme” from here on means the overall adaptive scheme procedure. 

(a) No Interpolation : For the linear problem, the behavior of the adaptive TVD schemes is 
similar to that of the classical shock-capturing methods. As opposed to the uniform grid case, 
the adaptive TVD schemes without interpolations behave rather poorly in terms of stability and 
allowable e values. See Figs. 1 1 and 1 2 of Yee & Sweby (1995b). The solutions refuse to settle 
down for larger At and/or smaller e. For the nonlinear problem, the behavior of the adaptive 
TVD schemes is similar to that of the uniform grid case. The range of allowable e and At in 
terms of stability and convergence rate and settling of the grid distribution are far better than in 
the linear problem. It appears that for problems with shocks, adaptive TVD schemes prefer no 
interpolation after each regridding. See Figs. 13 and 14 of Yee & Sweby (1995b). 

(b) With Interpolations : For the linear problem, as expected, both adaptive TVD schemes 
and adaptive classical schemes behave in a similar manner in terms of stability and convergence. 
Adaptive TVD schemes are less stable and have a smaller allowable range of e than the uniform 
grid case. Overall, adaptive TVD schemes behave far better than their counterparts without 
interpolation for the linear problem as can be seen in Figs. 1 1 and 12 of Yee & Sweby (1995b). 
For the nonlinear problem, the adaptive TVD schemes with interpolations behave like the 
classical shock-capturing method. They experience nonconvergence of the solution, and the 
grid distribution cannot settle down for a certain range of At. This can be seen in Figs. 13 and 
14 of Yee & Sweby (1995b). 

It is surprising to see the opposite behavior of the adaptive implicit TVD schemes for the two 
model PDEs with and without interpolation combinations, especially when the same physical 
parameters, discretized parameters and initial data were used. 


4.3. Mismatch in Implicit Schemes for Time-Marching Approaches 

When implicit methods are used to time-march to the steady states, it is a common practice 
in CFD to use a linearized and/or a simplified implicit operator (or mismatched implicit/explicit 
operators) to reduce operation count. The simplified implicit operators usually are first-order 
and the explicit operator retains higher order spatial accuracy. The simplified implicit operator 
might not be consistent with the original implicit scheme. It might also be nonconservative 
even though the original and/or the explicit operator are conservative. One popular formulation 
with these mismatched implicit/explicit operators is the “delta formulation” of Beam & 
Warming (1978) in conjunction with implicit LMMs time discretizations. The logic is that 
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if the solution converges, the explicit operator dictates the final accuracy of the steady-state 
numerical solution. As discussed in Section 2.4, these mismatched implicit/explicit operators 
might induce unwanted spurious dynamics into and/or reduce the convergence rate of the 
overall solution procedure. To illustrate just this point, we summarize some old work of Mulder 
& van Leer (1984) and the first author’s experiences (Yee 1986, 1989, 1990) in selecting the 
more desirable mismatch operators. These works use the delta formulation with a variant of 
the first-order upwind spatial discretizations for the implicit operator. The time discretization is 
the implicit Euler with the noniterative linearized form as discussed in Section 3.5.2. After the 
linearization in time and the drop in spatial discretization to first order, there are many ways to 
approximately evaluate the Jacobian matrix associated with the linearization for high-resolution 
and higher-order upwind shock-capturing schemes. Due to a space limitation, the readers are 
referred to the original papers for details or Yee ( 1989) for a summary. 

Mulder and van Leer studied two first-order upwind spatial discretizations (van Leer’s 
differentiable and Roe’s nondifferential forms) for the implicit operator and a first-order or 
second-order upwind spatial discretization for the explicit operator. They concluded that the 
differentiable first-order upwind implicit operator gives quadratic convergence for flow of an 
isothermal gas along an almost circular path through the stellar gravitational field of a rotating 
two-armed spiral galaxy. The governing equations are a system of hyperbolic conservation 
laws. With Roe’s nondifferentiable split flux-differences the iterations may get trapped in a 
limit-cycle. For a comparison of the two implicit operators, they used the same grid, physical 
parameters, explicit operator, initial data and numerical boundary treatment. 

Yee (1986) constructed conservative and nonconservative linearized forms of implicit TVD 
schemes. A rather detailed study on the convergence properties of these forms was performed. 
Using the same notation as the original paper, both the linearized nonconservative implicit (LNI) 
and the linearized conservative implicit (LCI) forms can be of first or second-order accurate 
in space. Numerical experiments in Yee et al. (1985), Yee (1986), Yee & Harten (1987) and 
Yee (1989) showed that both the second-order LNI and LCI forms are very unstable even if 
a very small time step is used. The first-order LNI and LCI perform far better. However, the 
first-order LCI form is more stable and has a better convergence rate than the LNI counterpart. 
The first-order LNI form has a higher chance of converging to a nonphysical solution. Also 
the residual sometimes stagnates or gets trapped in a limit cycle. These conclusions are based 
on comparing the two forms for a variety of one-dimensional and two-dimensional practical 
problems containing complex shock waves. For both LCI and LNI forms, the more compressive 
limiters, such as the superbee, are very unstable. The residual stagnates even with very small 
time steps. See Yee et al. (1990) for their performance for hypersonic computations. Grid 
refinement in this case does not improve the situation. Again, in order to isolate the cause of 
the convergence problem, the same grid, physical parameters, explicit operator, initial data and 
numerical boundary condition treatment were used. In passing, both the first-order LCI and 
LNI are heavily used in the CFD community with TVD schemes other than the Harten & Yee 
type (Yee 1989). This includes but is not limited to the various UNO, ENO and high-order 
upwind schemes for the explicit operator. 
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The above convergence properties using pseudo time marching approaches in conjunction 
with high- resolution implicit TVD schemes contributed to yet another kind of spurious numerics, 
which are very different from that in Sections 4.1 and 4.2. 


V. Spurious Dynamics in Unsteady Computations 

In Section 3.8, we cited some of the spurious numerics in transient computations. We have 
also learned that numerics can introduce and suppress chaos and can also introduce chaotic 
transients, the danger of relying on numerical tests for the onset of turbulence and chaos 
is evident. Here, some examples from CFD computations that exhibit analogous spurious 
behavior are illustrated. Sections 5. 1 , 5.2 and 5.4 are written by the authors of the original work. 
Laurence Keefe of Nielsen Engineering summarizes his unpublished work on chaotic transient 
computation that he performed in the late 1980’s. Shi Jin of the Georgia Institute of Technology 
summarizes his recent work on oscillations induced by numerical viscosities. Shi’s work has 
an important implication in spurious dynamics using time-marching approaches as well as 
slowly moving shock waves in transient computations. Section 5.3. summarizes the work of 
Brown & Minion (1995). It is concerned with spurious vortices in two-dimensional thin shear 
layer incompressible flow simulations using high-resolution shock-capturing methods. Bjorn 
Sjogreen of the Uppsala University, Sweden summarizes his work on convergence rate for 
systems of hyperbolic conservation laws in transient computations in Section 5.4. For additional 
results concerning other issues, see Moore et al. (1990), Corless (1992, 1994), Poliashenko & 
Aidun (1995) and Read & Thomas (1995). 

A brief background on chaotic transients and the significance and implications of Keefe’s 
work on numerical uncertainties are needed before presenting his results. Loosely speaking, 
a chaotic transient behaves like a chaotic solution (Grebogi et al. 1983). A chaotic transient 
can occur in a continuum or a discrete dynamical system. The concern here is ‘ ‘numerically 
induced chaotic transients” or “spurious chaotic transients”. One of the major characteristics 
of a numerically induced chaotic transient is that if one does not integrate the discretized 
equations long enough, the numerical solution has all the characteristics of a chaotic solution. 
The required number of integration steps might be extremely large before the numerical solution 
can get out of the chaotic transient mode. In addition, even without a numerically induced 
chaotic transient behavior, standard numerical methods usually experience a drastic reduction 
in step size and convergence rate near a bifurcation point (in this case the transition point) 
in addition to the bifurcation points due to the discretized parameters. See Section III for 
a discussion. Consequently, the possible numerically induced chaotic transient is especially 
worrisome in direct numerical simulations that are near the onset of turbulence. Away from 
the transition point, this type of numerical simulation is already very CPU intensive and the 
convergence rate is usually rather slow. Due to the limited computer resources, the numerical 
simulation can result in chaotic transients indistinguishable from sustained turbulence, yielding 
a spurious picture of the flow for a given Reynolds number. Thus, the danger of relying on 
numerical tests for the onset of turbulence is evident. 
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5.1. Chaotic Transients Near the Onset of Turbulence 
in Direct Numerical Simulations of Channel Flow 

Numerical simulations of wall-bounded turbulent shear flows have proven to be a powerful 
tool for investigating both the physics and mathematics underpinning these practically important 
flows. Boundary layers and turbulence are inescapable in aeronautical applications, and simple 
flows over a flat plate or within a channel provide an accessible arena in which to understand 
turbulent dynamics and to test ideas for its modification and control that can benefit the 
performance of real aircraft. 

Direct numerical simulations (DNS) of turbulent channel flow are the best developed of 
these techniques, with a better than 20 year history (Deardorff 1970, Schumann 1973, Moin & 
Kim 1982), and a relatively quick maturity (Kim et al. 1987) due to the favorable mapping of 
high accuracy spectral methods onto the geometry and known phenomenology of this flow. The 
physical situation is depicted in Fig. 5.1, where a flow is confined between planes at y = ±1 
and is driven in the ^-direction by a mean pressure gradient dp/dx. The flow is characterized 
by a Reynolds number Re = U^L/v, where Uoo is the mean centerline velocity, L is the 
channel half-height, and v is the kinematic viscosity. Within the channel the flow satisfies 
the incompressible Navier-Stokes equations and the no-slip boundary conditions are applied 
at the walls. In the particular calculations shown here these equations have been manipulated 
into velocity-vorticity form, where one integrates equations for the wall-normal velocity v, and 
normal vorticity tj and recovers the other two velocity components from the incompressibility 
condition and definition of tj. 
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Here the Hi contain the nonlinear terms in the primitive form of the Navier-Stokes equations 
and the mean pressure gradient. 
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Two experimentally observed facts figure strongly in the choice of numerical method used 
to integrate the Navier-Stokes equations in this geometry: the velocity increases extremely 
rapidly normal to the wall, and turbulent channel flows are essentially homogeneous in planes 
parallel to the wall. The first requires a concentration of grid points near the wall, and the 
second suggests use of a doubly periodic domain in planes parallel to the wall. Happily, a high 
accuracy spectral representation of the velocity field (u, v, tr) meets these needs: 

“ = £ £ £ (5.2) 

l m n 

where the Tj(y) are Chebyshev polynomials. The numerical problem then becomes dependent 
on a and 0 in addition to Re. For time advance mixed explicit-implicit methods are used. 
The nonlinear terms in the equations are advanced using second-order Adams-Bashforth or a 
low storage, third-order Runge-Kutta scheme (Spalart et al. 1991), while the viscous terms are 
advanced by Crank- Nicholson. The algorithm based on this spectral spatial representation and 
mixed time advance has been tested and used extensively, demonstrating an ability not only to 
reproduce experimental results, but to go beyond them in elucidating flow features not easily 
investigated in experiments. This code, and ones similar, (Handler et al. 1989, Jung et al. 1992) 
are currently being used to investigate a variety of turbulence control ideas suitable for turbulent 
drag reduction and (conceivably) separation control. 

One of the central problems in studies of wall bounded shear flows is the determination 
of when a steady laminar flow becomes unstable and transitions to turbulence. In dynamical 
systems’ terms, the Navier-Stokes equations always have a fixed point solution for low enough 
Reynolds numbers, but for each flow geometry the Reynolds number at which this fixed point 
bifurcates needs to be determined. In channel flow the fixed point solution (a parabolic velocity 
profile across the channel, u(y) = (1 — y 2 )) becomes linearly unstable at Re = 5, 772 (Orszag 
1971). However, since turbulence appears in experiments at much lower Reynolds numbers, 
it was conjectured that this bifurcation must be subcritical. Subsequent numerical solution 
of the nonlinear stability equations (Herbert 1976, Ehrenstein & Koch 1991) demonstrated 
this to be true, showing that limit cycle solutions with amplitude e branch back to lower 
Reynolds numbers before subsequently passing through a turning point and curving back 
toward higher Reynolds numbers. Thus for Reynolds numbers just above the turning point the 
flow equations have at least four solutions: the fixed point; two unstable limit cycles; and a 
chaotic solution (experimentally observed turbulence). Determining the location of the turning 
point in (a,0, e, Re) space is known as the minimum-critical-Reynolds-number problem, and 
its solution is by no means complete. 

One way to investigate the turning point problem is to perform DNS of channel flow for 
conditions believed to be near this critical condition. Beginning with a known turbulent initial 
condition from higher Reynolds number, one integrates the flow solver in time at the target 
Reynolds number to determine whether the flow decays back to the fixed point or sustains itself 
as turbulence. Although this may not be the most efficient way to bracket the turning point, 
it has the advantage that the peculiar dynamics of the flow near the turning point, whether in 
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decay or sustained turbulence, are observable, and this yields information about the path along 
which flows become turbulent at these low Reynolds numbers. 

Unfortunately the flow dynamics are very peculiar near the turning point, and extremely 
long chaotic transients are observed in the computations that make a fine determination of that 
point all but impossible by this method. This can be seen in Fig. 5.2, where a time history of 
the turbulent energy in a channel flow (energy above that in the laminar flow) is plotted for a 
Reynolds number of 2, 1 9 1 . To understand the time scale of the phenomenon some experimental 
facts need to be recalled. In typical experimental investigations of channel flow, the infinite 
transverse and streamwise extent of the ideal flow are approximated by studying flow in high 
aspect ratio (10-40) rectangular ducts that typically are 50-100 duct heights long. If times 
are non-dimensionalized by the centerline mean velocity ZToo and the duct half height L, then 
statistics on turbulence are gathered by averaging hot-wire data over intervals A tUoo/L ~ 200. 
In the simulations and figure the time scale is based on the friction velocity u r and L, where 
typically 15-20 u T ~ Uoo- Thus averaging over intervals A tu T /L ~ 10 should and does yield 
stable flow statistics that compare well with experiments. This can be seen in Figs. 5.3-5.5, 
where the near-wall velocity profile, cross- channel turbulence intensities, and Reynolds and 
shear stress distribution for the A tu r jL ~ 10 interval near the end of the transient, delineated 
by the arrows in Fig. 5.2, are shown. In each case they correspond well to available experimental 
data. Yet look at the time scale of the transient; it spans A tu r /L ~ 300, thirty times longer than 
the time needed to obtain stable statistics that would convince most experimentalists that they 
are viewing a fully developed turbulent channel flow. This is further complicated by the wide 
variation of the transient length, dependent upon both the grid resolution (number of modes in 
the spectral representation) and the linearly stable time step of the integration. In fact, for fixed 
(a, j3, Re) it is possible to obtain sustained turbulence for one time step, but see it rapidly decay 
to the laminar flow for another, lower value of the step. 

Extended chaotic transients near bifurcation points are not an unknown phenomenon; the 
"meta-chaos" of the Lorenz system is but one of many known examples. However, the 
practicalities of numerical computation in fluid dynamics usually interfere with ones ability 
to discern whether a transient, or sustained turbulence, is being calculated. The computations 
required to obtain the transient plot in Fig. 5.2 needed 40 hours of single processor time on 
a Cray XMP, some ten years ago. Such a small amount of expended time was only possible 
because the spatial resolution of the calculation was relatively coarse (32 x 33 x 32), in keeping 
with the large scales of the phenomena expected at these flow conditions. Higher resolution 
calculations (192 x 129 x 160) (Kim et al. 1987) at greater Reynolds numbers typically have 
taken hundreds of hours (~ 250) to barely obtain the A tu T /L = 10 averaging interval that is 
so inadequate for detecting transients. Because such calculations are so time consuming, one 
typically chooses an integration time step that is a substantial fraction of the linear stability 
limit of the algorithm, so as to maximize the calculated "flow time" for expended CPU time. 
However, it is clear from these transient results that this practice has some dangers when close 
to critical points of the underlying continuous dynamical system. Thus it appears that just as 
pseudo-time integration to obtain steady solutions can result in spurious results, genuine time 
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integration can result in chaotic transients indistinguishable from sustained turbulence, also 
yielding a spurious picture of the flow for a given Reynolds number. 

To conclude this section, here we give a feel for the number of integration steps required for 
the above numerical solution to get out of the chaotic transient mode. This transient calculation 
was performed using a time step of .0025. In these same units (time scaled by wall friction 
velocity u r and channel half height L) the transient calculation length was 409.6. Thus this 
calculation extended over 163,840 time steps. 

In Keefe et al. (1992) the dynamics of the computation in terms of “dimension” and 
“Lyapunov exponent” calculations were performed at a higher Reynolds number Re = 3200. 
For this higher Reynolds number, a transient occured using a smaller time step of 0.0015 and 
a transient calculation length of 644.52 or 429,680 time steps. To examine if chaos occured, 
Keefe et al. then took a central portion of that calculation at around 45,600 time steps and 
computed the Lyapunov exponent hierarchy from it. See Keefe et al. for additional details. 


5.2. Oscillations Induced by Numerical Viscosities in 1-D Euler Computations 


5.2.1. Introduction 

Earlier work has reported the difficulty of computing slowly moving shocks (Robert 1990, 
Woodward & Colella 1984), where first-order Godunov or Roe-type methods produce spurious 
long wave oscillations behind the shock and eventually ruin the downstream pattern. Here 
slowly moving means that the ratio of the shock speed to the maximum wave speed in the 
domain is much less than one. Several heuristic arguments, or improvements on the Riemann 
solver have been made in Arora & Roe ( 1 996), Jin & Liu (1996), Liu ( 1 995) and Woodward & 
Colella (1984). To investigate the dynamical behavior of shock-capturing methods for slowly 
moving shocks, we review the work of Jin and Liu (1996) using the traveling wave analysis 
and stability theory of discrete shocks. The goal is to carefully study this peculiar numerical 
phenomenon and to understand its formation and propagation, instead of solving this problem. 

Recall that the definition of a discrete traveling wave solution 4*", an approximation of 
U(xj, t n ),t„ = nAt, requires 

*r = (5.3) 

where 

*At/Ax=p/q, (5-4) 

with & the wave speed for some relative prime integers p and q. During a numerical calculation 
condition (5.4) may not hold. In other words, at different times the numerical viscous profiles 
correspond to different families of the traveling waves. It was shown by Jin and Liu that, 
at different time steps the numerical solutions correspond to different traveling wave profiles, 
unless (5.4) is met. The downstream oscillations are generated by such perturbations of the 
discrete traveling wave profile. The oscillations propagate along characteristics and behave 
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diffusively (decay in l* and Zoo)- The perturbing nature of the viscous shock (traveling wave) 
profile is the constant source for the generation of the downstream oscillations for all time. 

In their numerical experiments they also observed the periodic structure of the perturbing 
viscous shock profile. The period is essentially the time for the shock to propagate one spatial 
grid cell. 

In Section 5.2.2 a numerical example, using a Roe-type upwind scheme on a Riemann 
problem of the compressible Euler equations that admits slow shocks, is given. Among the 
numerical artifacts observed in this example are the momentum spikes at the shock, and the 
downstream oscillations. They are indeed numerical artifacts. The momentum spike is generated 
by the artificial numerical viscosity introduced in the continuity equation, which does not exist 
in the physical Navier-Stokes equations. The downstream oscillations are introduced by the 
dynamical behavior of the numerical viscous traveling wave profile. In Section 5.2.3 a traveling 
wave analysis on a “viscous isentropic Euler equations” formulation (Euler equations with a 
special linear viscosity term in both the continuity and momentum equations) is presented to 
show the existence of the momentum spike. This is compared with the momentum profile of the 
Navier-Stokes equations, which does not have the spike. In section 5.2.4 the dynamics of the 
downstream oscillations will be examined, and its relation with the stability and perturbation of 
the discrete shocks will be established. For details see Jin & Liu (1996). 


5.2.2. Numerical Solutions of a Slowly Moving Shock 


Consider the one-dimensional compressible Euler equations of gas dynamics, 


dtp + d m m = 0, (5.5a) 

dtm -|- d a (pu 2 + p) = 0, (5.5b) 

8 t E + d m [ti(Z7 + p)] = 0. (5.5c) 

Here p, u, m, p and E are respectively the density, velocity, momentum, pressure and total 
internal energy per unit volume with m = pu. For an ideal gas, the equation of state is given by 

P= (t - 1)(>® - ( 5 - 6 ) 

Let F be the flux function for (5.5) and let A denote the Jacobian matrix dF(U)/dU with 
U = ( p , m, E) T . The Euler equations (5.5) are hyperbolic with eigenvalues 

a 1 = u — c, a 2 = u, a* = u + c, (5-7) 


where c = y/jpfp is the local speed of sound. The right eigenvectors of A form the matrix 
R = ( R 1 , R 2 , R * ) given by 


R = 


1 

u — c 
H-uc 


1 1 
u u + c , 
\vl 2 H + ucj 


(5.8) 
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with H = -^—j/ + The inverse of R defines the left eigenvectors (£*, L 2 , L*) = R~ 1 of A 
by 

.«** + ?) l(-M-i) 

-R — I 1 - ii 6 ju —62 | , (5.9) 

!(*.-?) }(-*.«+!) |h 


with 


63 = 


7-1 


&, = &, 


u 


Let -Aj+ 1/2 be the Roe matrix satisfying (Roe 1981) 

Wi + .) - m,) = A j+l/1 (u j+I - 

By projecting Uj+i — Uj onto {J2 j+1 / 3 } one obtains the characteristic decomposition 


U 


3+1 




(5.10) 


(5.11) 


(5.12) 


J >=1 


In this decomposition the local characteristic variables a p +1 ^ 7 can be obtained using Roe’s 
average which perfectly resolves stationary discontinuities. 

Let Xj+ 1/2 be the grid points, Uj + i/ 3 be the pointwise value of U at * = Xj+ 1/2 . and Uj be 
the cell center value of U at Zj = (zj+ 1/2 + Zj-i/j)- A first-order upwind. Roe-type scheme 
(Smoller 1983), hereafter called UW, will be used to illustrate the phenomenon, although this 
phenomenon occurs for basically all shock-capturing methods. The numerical flux in UW is 
defined by 


3+1/2 


F(Vj) + F(U i+l ) 


;Sgn(A' + 1 /J )(if +1 - (5.13a) 


where 7 ^ is the component of F(Uj) in the p-th characteristic family. 


nv i ) = Y.+i R, i- ( 513b ) 

p = 1 


Example: A Riemann problem using the following initial data (Roberts 1990) will be solved. 
In this problem, there is a Mach-3 shock moving to the right with a shock speed s = 0.1096: 



3.86 


1 

U L = 

-3.1266 

27.0913 

if 0 < x < 0.5; Ur = 

-3.44 

8.4168 


if 0.5 < * < 1.(5. 14) 


Here 7 = 1.4 and the results at t = 0.95 are displayed in Fig. 5.6. The computation is carried 
out in the domain [ 0 , 1 ] with A* = 0 . 01 , A< = 0 . 001 . One can see that there is a momentum 
spike and the post-shock solution is oscillatory. 
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5.2.3. The Momentum Spikes 


A traveling wave analysis on the “viscous Euler equations” shows precisely the formation 
of the momentum spike. Consider the following “special viscous isentropic Euler equations” 
for density p and momentum m: 

d t p + d m m = ed mm p, (5.15a) 

d t m + d m 1 -p(p) = cd Bm m. (5.15b) 

. P 

Here the pressure p(p) = kp 7 for some constants k and 7 . This hyperbolic system has two 
distinct eigenvalues u ± c, where u = m/p is the velocity and c = y/^kp~<~ 1 is the sound 
speed. Although the true numerical viscosity is far more complicated than those appearing 
on the right-hand- side of (5.15), a study on (5.15) is sufficient for a full understanding of the 
numerical momentum spike. 

The traveling wave solution to (5.15) is examined. Let £ = where t is the shock speed. 


Then the traveling wave solution takes the form 

p(x,t) = <}>{(), m(®,<) = ^(£) (5.16) 

with asymptotic states 

^(± 00 ) = 4>±, V’Ci 00 ) — ^±- (5-17) 

The Rankine-Hugoniot jump conditions require 

-»(*+ -*-) + (*+-*-) = #, (5.18a) 

-«W+ “ + [x* 1 + jK^+) - ~r~ -K^-)l = *• (5.18b) 

Y<p+ <p- J 

First, the shock is assumed stationary (s = 0). It corresponds to the eigenvalue u — c. Then 
the jump condition (5.18) reduces to 

|±+J>W+)= j=- +*(*-), (5.19) 

and the Lax entropy condition gives 

0 < tt_ — c_ < u_ + c_, u + — c + < 0 < u+ + c + , (5.20) 

where u± — rj>± / <j>± and c = 


Applying the traveling wave solution (5.16) in (5.15) one gets the following ODEs: 

d(<f> = rf> — t/>_, (5.21a) 


th? tb* 

= — +p((t>) - -j p{4>-)- 


(5.21b) 
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This system has two fixed points: V+ = on the right and V_ = on the 

left in the phase plane of ( <f > , r/>). It has two distinct eigenvalues, Aj = u — c and Aj = u + c, 
with corresponding eigenvectors iZ, = (l,u - c) T and R 3 = (l,u + c) r . By the entropy 
condition (5.20), V + is a saddle point with a stable manifold on R x , and VL is a source. Thus a 
heteroclinic orbit O will connect VL and V+ in the direction of Ri (Smoller 1983), as shown 
in Fig. 5.7. (In Fig. 5.7 R± and Rf are the two eigenvectors at V± respectively). The orbit 
O is smooth, and d^<f> is not identically zero if <f>— ^ 4 >J r- Thus (5.21a) implies that t/> is not a 
constant. Moreover, whenever </>(£) connects <f>- and <f>+ with a monotone profile, d(<f> becomes 
a spike. Thus ip = d^<f> 4- xp- is a spike. 

For a nonstationary shock, the traveling wave solution (5.16) applied to (5.15a) gives 

d(4> = -s{4> - <p-) + ip - ip-, (5.22) 

or 

rj> = s<f> + dt<j>+ {ip- - »(}>-). (5.23) 

Hence ip is a superposition of a monotone profile »<p with a spike corresponding to d(<p. When s 
is small (for a stationary or a slowly moving shock), the monotone profile »<p becomes small and 
the spike term d^<p dominates. Thus the shock profile of ip is a non-monotone spike. Therefore 
the spike is usually generated in a stationary or slowly moving shock, as shown in the earlier 
examples. For a strong shock the monotone profile s<p dominates so the shock profile of the 
momentum is monotone. 

Since the more physical viscous shock profile is determined by that of the Navier-Stokes 
equations, the viscous profile of the isentropic Navier-Stokes equations is now studied and 
compare it with that of the viscous Euler equations (5.15). The isentropic Navier-Stokes 
equations are 

d t p + d m m = 0, (5.24a) 

d t m + d m ?y+p{p) = e0„f™Y (5.24b) 

Applying the traveling wave solution (5.16) in (5.24) and again assuming the shock speed 
s = 0, one obtains the following ODEs: 

ip = (5.25a) 

(5.25b) 

tb 2 

Equation (5.25a) shows that ip is a constant and thus contains no spike. Let /(</>) = ~ p(<p), 

then (5.25b) becomes 

e ( (^) = m - M-). (5.26) 
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Since is a strict convex function and - f{4>~) is always negative 

between <f>- and <f>+. Therefore, 6^(1/ <f>) does not change sign by (5.26). This implies the 
monotonicity of l/<f> or <f>. 

When < 7 ^ 0, applying the traveling wave solution (5.16) in (5.24a) gives 

= s<f> + (^_ — s<f>-). (5-27) 

Thus whenever <)> is monotone so is t/>. This excludes the possibility of a momentum spike for a 
moving viscous shock in the Navier-Stokes equations. 

In conclusion, even if the viscous profile of <f> of the Navier-Stokes equations (5.24) could 
be similar to that of the viscous Euler equations (5.15), the profiles of xj> may be significantly 
different. Since the physically relevant solution of the Euler equations is considered to be 
the zero viscosity limit of the Navier-Stokes equations, the momentum spike appearing in the 
viscous Euler equations is totally nonphysical. 

By examining the interrelation between the viscous Euler and the Navier-Stokes equation 
one can come up with a change of variables that recovers the Navier-Stokes equations in 
an asymptotic sense from the Euler equations. This also motivates a numerical change of 
variable from the cell-center or cell average momentum to the mass flux, which eliminates the 
momentum spike exactly. For details see (Jin & Liu 1996). However, what is more catastrophic 
is the downstream oscillations, which cannot be easily eliminated, as will be studied next. 

5.2.4. The Downstream Oscillations 

Note that almost all shock-capturing methods are in conservative form. Due to the 
conservation of momentum, the total mass of momentum carried by the spike profile should be 
compensated by an equal amount of momentum mass elsewhere. This explains the formation 
of the downstream waves. 

Figure 5.8 shows the result of UW1 for the previous example after 5 time steps to illustrate 
the formation of the momentum spike and downstream wave. As the density is smeared, the 
momentum forms a spike and a downstream wave. The spike and the downstream wave carry 
the same mass so the total momentum is conserved. 

In order to demonstrate that the downstream oscillations propagate along characteristics and 
are diffusive, the Roe decomposition (5.12) is used, where a p represents the component of 
Uj +i — Uj in thep-th characteristic family. The numerical “characteristic” variable is defined 
as 

# = XX+./> A *- ( 5 - 2 «) 

»<j 

A distinction between the dispersive oscillations associated with center difference schemes and 
the downstream oscillations studied here is that the latter lie only in their own characteristic 
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family. For example a wave which appears in /?* does not appear in (3 q for p ± q. These can be 
seen in Fig. 5.9 where each wave moves away with the corresponding characteristic speed, and 
behaves diffusively (spread out and decay). 

Figure 5.10 displays the time evolution of the momentum profile of the first-order Roe’s 
scheme for the same example. One can see that the spike (viscous) profile keeps fluctuating in 
an 0 ( 1 ) manner, causing the downstream oscillations for all time. The diffusive nature of the 
downstream oscillations is evident in the picture. 

Figure 5.1 1 shows the peak of the momentum spike as a function of time. The fluctuating 
nature of the spike is evident. The more the mass of the spike profile varies the more strongly the 
diffusion waves emerge for momentum conservation. Interesting is that the peaks are periodic, 
with the duration of each period agreeing with the time for the shock to move one grid point. 
Consequently, the discrete shock profile is stable only modulo this period. However, within 
each period it is fluctuating, which becomes the source of the new downstream waves. 

Recall that the definition of a discrete traveling wave solution an approximation of 
U(xjyt n ),tn = nAt, requires = $j_ np , where sAt/Ax = p/q for some relative prime 
integers p and q. The stability of such a discrete shock for the Lax-Friedrichs scheme was 
established by Jennings (1974) for scalar equations and by Majda and Ralston (1979) and Liu 
& Xin (1993) for nonlinear systems. The periodicity of the momentum peaks in Fig. 5.1 1 
shows the stability of the discrete traveling wave solution for these schemes’ modulus the 
time for the shock to travel one grid point This is because, when a At « Ax, there exists a 
sufficiently large q such that |g(«A<) - Ase| < sAt, or \»AtjAx — l/q\ < 2/q 2 . However, 
within each period the numerical shock layer is unsteady and corresponds to different traveling 
wave profiles, which become the source of the new downstream waves in all time for these 
schemes. 

A scheme that completely eliminates the downstream oscillations in later time should have 
a steady viscous profile beyond the initial formation of the spike; i.e., the momentum spike 
peak should remain a constant in later time. However, this is impossible as long as the shock is 
moving and it takes many times for the shock to move to the next cell. 

The discrete shock profile perturbs even when the shock does 'not move slowly. Thus 
the downstream oscillations exist even for fast shocks. However, in the fast shock case the 
momentum profile is monotone and thus does not leave much room for the shock profile to 
perturb. In other words, each perturbation does not change the mass of the viscous profile 
much, and the downstream errors become negligible. For slow shocks the momentum profile 
has a spike, which increases the mass of the viscous profile and the relative mass change in each 
perturbation. Therefore, the downstream errors become more significant. This also illustrates 
why the downstream errors in the density are far less significant. Since the density is monotone, 
the relative change in the mass of the viscous profile is smaller than that of the momentum. 

In summary, although each family of the downstream waves decay time- asymptotically, the 
perturbing spike or viscous profile is a constant source for the generation of new downstream 
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waves, causing the downstream noise for all time. Higher order methods use higher order 
interpolations, which amplify the noise and exhibit rich but spurious post-shock structures. 

5.2.5. Discussions and Conclusions 

As studied in (Jin & Liu 1996), similar behavior occurs in schemes that are of monotone, 
TVD or ENO type. Note that all these monotonicity theories are established only for scalar 
equations, or linear systems via the characteristic variables. For nonlinear systems there are no 
global characteristic variables. Thus these methods are usually extended to nonlinear systems 
using the idea for linear systems; i.e., via the so-called local characteristic decomposition (using 
the Roe matrix example). Since there is no theory for the monotonicity of these methods for 
nonlinear systems, it is not surprising to see the non-monotone behavior represented by the 
spike and the downstream oscillations reported here. As pointed out by Jin & Liu (1996), 
to fully solve this problem, instead of applying scalarly monotone, TVD or ENO scheme to 
nonlinear systems, one may need a method that is systematically “monotone, TVD or ENO’’. 
One may also need to choose numerical viscosity properly so it mimics the physical viscosity 
of the Navier-Stokes equations. The ultimate goal is to have a scheme that not only provides a 
high resolution but, more importantly, has a more stable viscosity profile. 

The novelty of the work of Jin and Liu is that they are the first to use the traveling wave 
analysis to prove the non-monotonicity of the solution for nonlinear systems, and to link the 
downstream oscillations to the stability of discrete traveling wave profile. To really understand 
the behavior of shock-capturing methods for nonlinear systems, and to ultimately design 
nonoscillatory schemes for nonlinear systems, good theories for both inviscid and viscous 
nonlinear systems need to be developed. This remains an open and challenging research subject 
for the future. 

5.3. Spurious Vortices in Under-Resolved Incompressible Thin Shear Layer Flow Simulations 

Brown & Minion (1995) performed a thorough study of a Godunov-projection method and a 
fourth-order central difference method for the two-dimensional incompressible Navier-Stokes 
equations as a function of the resolution of the computational mesh with the rest of the physical 
and discretized parameters fixed. In the authors’ opinion, this is a good example of isolating 
the cause of the spurious behavior. The physical problem is a doubly periodic double shear 
layer. The shear layers are perturbed slightly at the initial time, which causes the shear layer to 
roll up in time into large vortical structures. For a chosen shear layer width that is considered 
to be thin and a fixed perturbation size, they compared the solution for four different grid sizes 
(64 x 64, 128 x 128, 256 x 256, 512 x 512) with a reference solution using a grid size of 
1024 x 1024. For the 256 x 256 grid, a spurious vortex was formed midway between the 
periodically repeating main vortex on each shear layer. The 128 x 128 solution showed three 
spurious vortices along the shear layer. The spurious vortex disappeared with a 512 x 512 
mesh. They also disabled the limiters (a strictly upwind Fromm’s method), and found the 
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behavior to be similar. They also cited other work using Lax-Wendroff type and ENO schemes 
where similar behavior was observed. They concluded that the spurious vortex is the artifact of 
underresolution of the grid. Linking this behavior with a re-interpretation of their conclusion 
using nonlinear dynamics, we interpret their observation as follows. For the particular grid size 
and time step combination, stable spurious equilibrium points were introduced by the numerics 
into a portion of the flow field while the major portion of the flow field was predicted correctly. 
In other words, the spurious vortices are the solution of the discretized counterpart for that 
particular range of grid size and time step. The number of stable spurious vortices is a function 
of the grid size. As the grid spacing decreases, the spurious equilibriums gradually become 
unstable and the numerical solution mimics the true solution. 


5.4. Convergence rate for systems of hyperbolic conservation laws 

In Section 4.2.1 the system (4.7a), (4.7b) was solved with the initial data 

f 1 * < 0 

txo(z) = | ft » tfot*) = 1* (5.29) 

l — 1 * > 0 

The exact solution is then independent of time. In order to see the dynamical behavior of the 
error propagating from the shock, here a time dependent problem is solved instead. This can be 
accomplished by changing to a different initial data. Sjogreen (1996) use the initial data 

f 1 * < 0 

tto(x) = < i g »o(x) = sin** (5.30) 

for the same system of partial differential equations (4.7). With these data, the exact solution is 

u(x,f) = uo(x) v(x 7 1) = sinir(x — <). (5.31) 

He solves the problem numerically on the interval — 1 < * < 1, using 100 uniformly 
distributed grid points. The numerical method is again a fully second order accurate scheme of 
ENO type in space, and a second order accurate Runge-Kutta method in time. The exact value 
v(— 1,<) = sinjr(— 1 — t ) is given on the left boundary and second order accurate extrapolation 
is used at the boundary x = 1. 

Figure 5.12 shows the u and v -components of the solution at the time 0.265958, computed 
on 100 grid points. There is an error in the v component, similar to what can be seen in Figure 
4.2. However here it is harder to see it in the plot, since the sine wave solution forces a different 
scaling of the y-axis. 

A better way of seeing the error is in Fig. 5.13, where the logarithm of the pointwise error in 
the v component is plotted against x. The three plots correspond to 30, 100, and 200 time steps 
respectively. A domain of first order grid convergence is propagating out from the shock, seen 
as a flat region of large error in the plots. The first order convergence in this region was, in this 
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case, obtained from numerical experiments, although analytical examples are also possible (see 
Engquist & Sjogreen (1995)). 


VI. Concluding Remarks 

The need for the study of dynamics of numerics is prompted by the fact that the type 
of problem studied using CFD has changed dramatically over the past decade. CFD is also 
undergoing an important transition and it is increasingly used in nontraditional areas. But even 
within its field, many algorithms widely used in practical CFD applications were originally 
designed for much simpler problems, such as perfect or ideal gas flows. As can be seen in 
the literature, a straightforward application of these numerical methods to high speed flows, 
nonequilibrium flows, advanced turbulence modeling or combustion related problems can lead 
to wrong results, slow convergence, or even nonconvergent solutions. The need for new 
algorithms and/or modification and improvement to existing numerical methods in order to 
deal with emerging disciplines is evident. We believe the nonlinear dynamical approach for 
CFD can contribute to the success of this goal. The first step toward achieving this goal is 
to understand the nonlinear behavior, limits and barriers, and to isolate spurious behavior of 
existing numerical schemes. 

We have revealed some of the causes of spurious phenomena due to the numerics in an 
attempt to improve the understanding of the effects of numerical uncertainties in CFD. We 
illustrated with practical CFD examples that exhibit similar properties and qualitative behavior 
as elementary examples in which the full dynamical behavior of the numerics can be analyzed. 
We have also shown that guidelines developed using linearization methods are not always valid 
for nonlinear problems. Even well-informed use of conventional methods may lead to nonsense 
on unconventional problems. We have gained an improved understanding of long time behavior 
of nonlinear problems and nonlinear stability, convergence and reliability of time-marching 
approaches. We have learned that numerics can introduce and suppress chaos and can also 
introduce chaotic transients and the danger of relying on numerical tests (e.g., direct numerical 
simulation) for the onset of turbulence and chaos is evident. The nonlinear phenomenon and 
spurious behavior exhibited by the numerics in solving genuinely nonlinear problems reveal 
many of the limitations, challenges and barriers in CFD. We believe the knowledge gained 
so far has already provided some improved guidelines for overcoming the spurious behaviors 
without resorting entirely to the tuning of computational parameters. Since standard procedures 
such as using physical guidelines, grid refinement, halving of the time step, and using more than 
one scheme to assure the quality, reliability and the integrity of the numerical solution for stiff 
and strongly nonlinear problems are not foolproof and/or not always possible, the dynamical 
systems approach can be a viable complement. Before additional theories are established, we 
conclude that the safest route is to have some understanding of the dynamical behavior, limits 
and barriers of the numerical method being used. 

As can be seen, recent advances in dynamics of numerics showed the advantage of adaptive 
step size error control for long time integrations of nonlinear ODEs. Although much research 
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is needed to construct suitable yet practical similar adaptive methods for PDEs, these early 
developments lead our way to future theories. There remains the challenge of constructing 
adaptive step size control methods that are suitable yet practical for time marching to the 
steady states for aeronautical CFD applications. Another even more challenging area is the 
quest for an adaptive numerical scheme that leads to guaranteed and rapid convergence to the 
correct steady-state numerical solutions. These two key challenges are particularly important 
for CFD. We conclude the paper with the following guidelines to minimize spurious dynamics 
in time-marching to the steady state. 

Some Guideline » to Minimize Spuriout Dynamic < : Due to the spurious dynamics intro- 
duced by the numerics, one usually will not be able to map out the complete numerical basins 
of attraction and bifurcation diagrams for the entire problem in practical situations. Only in 
isolated situation with a particular physical problem and numerical method combination such 
as the one studied in Lafon & Yee (1991, 1992) are continuation methods able to locate all of 
the essential spurious branches of the bifurcation curves. 

On the other hand, continuation methods are widely used in dynamical systems when one 
wants to understand certain properties of key branches of the bifurcation curve, especially if 
one knows (or can ascertain by other means) a starting solution on that particular branch. For 
example, in the Taylor-Couette flow problem, extensive use is made of continuation methods to 
map out the critical Reynolds number when the flow behavior undergoes drastic changes in flow 
patterns, since in this case we know how the flow behaves for the low Reynolds number case. 
Another example is in Bailey and Beam (1991) who used it to study the hysteresis behavior of 
the flow of an airfoil in terms of angles of attack for the steady PDEs. In this case, the flow 
behavior is readily obtained for low angles of attack (before hysteresis). Most of the use of the 
continuation method so far is focused on elliptic PDEs or steady PDEs. These studies seldom 
address the possibilities of spurious dynamics due to the numerics, especially for IB VPs using 
time-marching approaches. It is remarked that a shortcoming in association with solving the 
steady PDEs is that a small radius of convergence or nonconvergence of the numerical solution 
is often encountered even with the aid of multigrid, preconditioners, and/or relaxation methods, 
especially when the PDEs are of the mixed type (e.g., the steady inviscid supersonic flow over 
a blunt body). In addition, the solutions obtained do not distinguish whether the steady states 
are stable or not. 

Besides the study in Lafon & Yee (1991, 1992), here we propose a further step of applying 
this technique to the discretized counterparts of the time-dependent PDEs in order to avoid 
spurious asymptotes due to unknown initial data. The idea relies on the knowledge of a known 
or a reliable numerical solution on the correct (non-spurious) branch of the bifurcation curve 
as a function of the physical parameter of interest. The logic is that if one starts on the correct 
branch, one avoids getting trapped on any of the spurious branches. Also the issue of unknown 
initial data related to time-marching approaches is avoided or can be minimized. Details of the 
approach and numerical examples will be reported in a forthcoming paper. Here we will give a 
short narrative summary of the procedure. 
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In many fluid problems where the solution behavior is well known for certain values of the 
physical parameters, but unknown for other values. For these other values of the parameters, 
the problem might become very stiff and/or strongly nonlinear, making the available numerical 
schemes (or the scheme in use) intractable. In this situation, continuation methods in bifurcation 
theory can become very useful. If possible, one should start with the physical par am eter of a 
known or reliable steady state (e.g., flow behavior is usually known for low angles of attack 
but not for high angles of attack). One can then use a continuation method such as the pseudo 
arclength continuation method of Keller (1977) (or the recent development in this area) to 
solve for the bifurcation curve as a function of the physical parameter. The equations used are 
the discretized counterpart of the steady PDEs or the time-dependent PDEs. If time-marching 
approaches are used, a reliable steady-state numerical solution (as a starting value on the correct 
branch of the bifurcation curve for a particular value of the physical parameter) is assumed. 
This starting steady-state numerical solution is assumed to have the proper time step and initial 
data combination and to have the grid spacing fine enough to resolve the flow feature. The 
continuation method will produce a continuous spectrum of the numerical solutions as the 
underlying physical parameter is varied until it arrives at a critical value p c such that it either 
experiences a bifurcation point or fails to converge. Since we started on the correct branch of 
the bifurcation curve, the solution obtained before that p c should be more reliable than if one 
starts with the physical parameter in question and tries to stretch the limitation of the scheme. 
Note that by starting a reliable solution on the correct branch of the bifurcation curve, the 
dependence of the numerical solution on the initial data associated with time-marching methods 
can be avoided. 

Finally, when one is not sure of the numerical solution, the continuation method can be used 
to double check the numerical solution. This approach can even reveal the true limitations of 
the existing scheme. In other words, the approach can reveal the critical physical parameter 
for which the numerical method breaks down. On the other hand, if one wants to find out the 
largest possible time step that one can use for a particular problem and physical parameter, one 
can also use continuation methods to trace out the bifurcation curve as a function of the time 
step. In this case, one can start with a small time step with the correct steady state and observe 
the critical time step as it undergoes instability or bifurcation. 
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Scheme 

Fixed points 

Stable range 

Explicit Euler 

1 

0 <r<2 

Modified Euler 

1 

0<r<2 


1 + 2/r 

0<r<-l + V5« 1-236 


2/r 

2 < r < 1 + V5 = 3-236 

Improved Euler 

1 

0 < r <2 


[2 + r±Vr 2 -4] 
2 r 

2<r<V8 = 2-828 

Heun 

1 

0< r < 1 + (Vl7 + 4) ,/3 — (Vl7 + 4) _l/3 = 2-513 


♦ 

4-9137 <r <4-9552 


* 

6-4799<r< 6-4853 


* 

6-74405 <r< 6-74575 

R-K 4 

1 

0<r<? + (W + 5V29)‘ /3 

+ (W-5V'29) 1,, = 2-785 


♦ 

2-785 <r <3-4156 


* 

2-746 <r <3-456 


Table 3.1. Fixed points of Runge-Kutta methods for du/dt = atx(l - «). 


Scheme 

Fixed points 

Stable range 

Explicit Euler 

\ 

0<r<8 

Modified Euler 

1 

2 

0< r <8 


i[l + Vl — 32/r] 

32 <r <32-014067 


i(l±Vl-8/r] 

8 < r < 4(1 + V5) = 10-928 


«[3 — Vl — 32/r] 

32 <r< 32-014067 

Improved Euler 

1 

2 

0<r<8 


1 ± Vl - 8 It 
2 

8 < r < 12 


1 Vl-12/r Vl+4/r 

2 4 1 4 

12 < r < 4(1 + V6) * 13-798 


1 Vl - 12/r Vl + 4/r 
2 + 4 ± 4 

12 < r < 4(1 + V6) * 13-798 

Heun 

5 

0< r <4(1 + (Vl7 + 4) ,/3 



-(Vl7 + 4)- ,/3 )= 10-051 

R-K 4 

1 

2 

0 < r < + M/g + |V29) ,/3 



+ 4(^-^V29) ,/3 « 11-14 


Table 3.2. Fixed points of Runge-Kutta methods for du/dt = ati(l - u)(0.5 - u). 
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S(u) 

(3.11) 

(3.17) 

(3.18) 

(3.19) 

(3.21) 

“(l -«) 

0 

2 

2 

6 

14 

“(2 - w)(l -«) 

0 

6 

6 

24 

78 


Table 3.3. The number of spurious fixed points of Runge-Kutta methods 
for (3.7a) and (3.23). 


Equation 

Period 2 orbits 

Stable range 

u' = au{ 1 - u) 

r + 2 ± Vr 2 — 4 

2r 

2<r<V6 = 2-4495 

u' = au(l - u)(i ~ u) 

K, 

00 

1 

-H 

8 < r < 12 


1 Vl-12/r Vl+4/r 

2 4 ± 4 

12 < /■ < 14 


1 Vl-12/r Vl+4/r 

2 + 4 ± 4 

12 < /• < 14 


Table 3.4. Period 2 fixed points of the explicit Euler method. 
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scheme 

Ci 

c i 

E-0 

0.65 

0.9 

minmod 

0.5 

0.75 

superbee 

— 

0.7 

van Leer 

0.6 

0.7 

van Albada 

0.6 

0.7 


Table 4.1. Convergence regions for the explicit schemes. 


scheme 

converges to x\ 

converges to x 2 

E-0 

c < 11 

c > 22.5 

minmod 

c < 11 

c > 21.5 

superbee 

c < 11 

— 

van Leer 

— 

— 

van Albada 

c< 11 

c > 22.5 


Table 4.2. Convergence regions for the implicit schemes. 



at stable shock 

at unstable shock 

scheme 

HI 

Y 

stable 

X 

Y 

stable 

explicit E-0 

0.36 

-0.47 

At < 0.057 

mm 

mm 






mm 







Ell 



implicit E-0 

0.37 

-0.48 

VA* 

■ssa 

BE9 

At > 1.0585 





WSM 

fli 

VA< 





-0.29 

-0.53 

VAt 


Table 4.3. Analytical fixed points of the reduced system. 
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7 • 


IS 


J T n r 

7 a • 
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Fig. 3.3. Stable fixed points of periods 1, 2, 4, 8 of die modified Euler method for the ODE 
du/dt = au(l — u)(6 — u). 
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r i l l l I 


1.8 2.0 2.2 2.4 2.6 2.8 3.0 

r a* aAt 


DIAGRAM LOOKS THE SAME REGARDLESS OF I.C. FOR ALL u°>0 
CHAOS WINDOWS NEAR: 2.627, 2.634, 2.738, 2.828, BELOW 3 


Fig. 3.4. Bifurcation diagram of the explicit Euler method for the logistic ODE (a > 0). 
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Fig. 3.6. Bifurcation diagrams of the modified Euler method for the logistic ODE (a > 0). 
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Bifurcation Diagrams & Basins of Attraction 

u' = au (1-u) 


Modified Euler 


Improved Euler 




Kutta 



R-K 4 




Fig. 3.9. 
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Bifurcation Diagrams & Basins of Attraction 

u’ = au (1-u) 
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Bifurcation Diagrams & Basins of Attraction 

u’ = au (1 -u) 
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Fig. 3.11. 







Bifurcation Diagrams & Basins of Attraction 

u’ = a u (1-u) 


O 
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Fig. 3.13. 




Bifurcation Diagrams & Basins of Attraction 

u’ = au(1-u) 
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Fig. 3.14. 




Bifurcation Diagrams & Basins of Attraction 

u’ = a u (1-u) 
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Fig. 3.15. 






Bifurcation Diagrams & Basins of Attraction 

u’ = au(1-u) 
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Fig. 3 16. 




Bifurcation Diagrams & Basins of Attraction 
u’ = a u (1 - u) (0.2 - u) 



ill! 




Bifurcation Diagrams & Basins of Attraction 
u’ = a u (1 - u) (0.2 - u) 



mi . mi 




Bifurcation Diagrams & Basins of Attraction 
u’ = a u (1 - u) (0.2 - u) 
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Fig. 3.19. 




L2 norm of U Spectral radius 




Rt 3.20. Comparison for discretisation UF1UI/EE and „ = 7 of (a) linearized stability analysis 
(around u* = 2) and (b) nonlinear solution behavior with initial data U 
(• : steady-state solution and other asymptotic solution, blank space : verger. 

solution). 
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L2 norm of U Spectral radius 




P2 = aAt 


Fig 3 21 Comparison for discretization UPWI/ME and pi = 7 of (a) Lnoarmad stabity anal- 
' vsis (around u* = 2) and (b) nonlinear solntion behavior with initial data V (• 
steady-state solntion and other asymptotic solution, blank space . divergent solution). 





0 12 3 4 

PI 

Fig. 3.22. Average wave speed W versus pi of the numerical solution with explicit Euler time 
discretization and spatial discretization UP1PW. 
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CHA/EE CHA/ME 



Fig. 3.23. Numerical basin of attraction for discretizations (a) CHA/EE, (b) CHA/ME, (c) 
CHA/LIE and (d) CHA/LT with J = 4, ft = 0.1 andft = 0.5 (A : exact steady-state 
solution, ■ : spurious steady-state solution, • : other spurious asymptotic solution, 
blank space: divergent solution). 
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CHA/LIE CHA/LT 



Fig. 3.24. Numerical basin of attraction for discretizations (a) CHA/LIE and (b) CHA/LT with 
J = 4, pi = 0.1 and p 2 = 3 (CFL = 0.3) (A : exact steady-Btate solution, ■ : spurious 
steady-state solution, • : other spurious asymptotic solution, blank space: divergent 
solution). 


CHA/LIE CHA/LT 



Fiu. 3.25. Numerical basin of attraction for discretizations (a) CHA/LIE and (b) CHA/LT with 
J ~ 4, pi = 0.1 and pj = 10 (CFL = 1) (A : exact steady-state solution, ■ : spurious 
steady-state solution, • : other spurious asymptotic solution, blank space: divergent 
solution). 
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Fig. 4.4. Error in entropy (left) and convergence exponent (right) of (4.7) for a second-order 
UNO scheme. 


Ull 


fixed 



Fig. 4.5. Grid points of the reduced Embid et al. problem. 


3 ' 
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slocity Profile 



les. o; "corrected" data of Eckelmann (1974); 

- : "law of the wall" : u + = y + , u + = 2.5 ln(y + ^ + 5.5 


-nee Intensities 



vA 5 


jations normalized by wall shear velocity, 
□ss-channel coordinate normalized by 5 = 2L. 
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Reynolds Shear Stress 



Fig. 5.5. Reynolds shear stress normalized by wall shear velocity. , -liv; 

-uv + (1 /Re)du/dy, , Total shear stress for fully developed channel 


Density 


Momentum 



Fig. 5.6. A slowly moving shock at t = 0.95 computed by UW1 using Az = 0.01 and 
A t = 0.001. 
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<t>+ <p 

g. 5.7. A sketch of the phase portrait of Eqs.(5.21). 


Momentum 




Fig. 5.9. The downstream waves in “characteristic” variables. Note that each wave belongs 
only to one characteristic family and diffuses. 





Fig. 5.10. Time evolution of the momentum spike and the downstream waves by the Roe s 
lst-order upwind scheme for t £ [0.5, 0.8]. For better visualization these graphs are displayed 
upside down. The diffusive nature of the downstream pattern is apparent. 



Spike Peak 



Time 


Fig. 5.11. Time evolution of the peak of the momentum spike by the Roe’s lst-order upwind 
scheme 
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